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ABSTRACT 


In  this  thesis  the  optimal  torque  and  voltage  control  of  a 
turbo-generator  connected  to  an  infinite  bus  is  considered. 

The  generator  is  controlled  through  a  linear  feedback  of  the 
state  variables.  The  feedback  parameters  are  obtained  by  solving 
a  two  point  nonlinear  boundary  value  problem.  The  values  to  be 
obtained  for  these  parameters  depend  on  the  strength  and  duration 
of  the  disturbance  since  the  model  is  nonlinear  contrary  to  the 
usual  feedback  control  of  a  linear  model. 

The  methods  of  solution  are  to  cast  the  system  into  the 
Pontryag  in  ’  s  minimum  principle  and  the  gradient  descent  method. 
The  model  used  includes  the  transfer  functions  of  the  governor, 
the  turbine  and  the  exciter. 
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CHAPTER  I 


INTRODUCTION 

1 . 1  Preliminary  Remarks 

Modern  trends  in  power  systems  towards  the  use  of  generating 
units  of  larger  capacity  and  lower  specific  inertia,  and  longer 
transmission  lines  at  higher  voltages  tend  to  decrease  the  transient 
stability  and  steady  state  stability  limits  of  the  system.  The 
problem  of  improving  the  steady  state  stability  has  been  studied 
and  the  various  methods  proposed  have  been  shown  to  provide  an 
adequate  margin  of  stability  in  the  leading  power  factor  region  of 
operation  even  under  low  loads. 

As  the  steady  state  stability  limits  are  constantly  improved, 
the  operation  of  the  power  system  tends  to  become  increasingly 
constrained  by  its  transient  stability  limits.  Methods  which  have 
been  proposed  for  increasing  the  transient  stability  margin  act 
generally  by  controlling  the  output  of  the  generator  during  transient 
cond it ions . 

Regulation  of  the  power  input  to  the  prime-mover  will  undoubtedly 
produce  a  direct  effect  on  the  rotor  swings  of  the  generator  under 
transient  conditions  and  can  be  achieved  in  the  case  of  a  thermal 
set  by  control  of  the  steam  valve.  However,  due  to  the  presence  of 
entrained  steam  storage  effects  in  the  various  parts  of  the  turbine, 
in  the  reheater  and  in  the  connecting  pipes,  the  maximum  benefit 
from  this  type  of  control  are  not  generally  achieved.  Moreover, 
the  centrifugal  type  of  governor  and  control  equipment  are  designed 
mainly  to  prevent  over speeding  of  the  turbine  under  conditions  of 
major  load  rejection  or  line  switching. 
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From  these  observations,  to  obtain  the  benefits  of  turbine  regulation 
during  the  transient  period,  it  is  necessary  to  consider  its  effect 
simultaneously  with  that  of  the  automatic  voltage  regulator. 

The  method  of  approach  is  through  the  use  of  dynamic  optimization 
techniques  to  determine  an  optimal  feedback  control  policy  subject  to  a 
prescribed  performance  criterion.  For  this  purpose  either  the 
calculus  of  variations  or  the  Pontryagin ’ s  minimum  principle  could  be 
used.  The  criterion  for  optimum  transient  performance  is  taken  as  an 
integral  quadratic  index  which  includes  soft  constraints  on  the  control 
variables . 

1.2  Scope  of  the  Thesis 

In  this  thesis  a  general  nonlinear  model  of  a  turbo-alternator 
connected  to  an  infinite  bus  is  considered.  Control  of  the  turbo¬ 
alternator  is  effected  through  control  of  field  voltage  and  turbine 
torque.  In  addition,  it  is  proposed  to  include  the  transfer 
functions  of  the  governer,  the  turbine  and  the  voltage  regulator. 

Inclusion  of  these  transfer  functions  provides  a  more  realistic 
model  of  the  machine  and  its  control  system,  however,  the 
dimensions  of  the  system  are  increased. 

The  machine  is  disturbed  from  the  steady  state  by  a  torque 
pulse  applied  directly  to  the  shaft  andwhile  the  alternator  is 
controlled  through  a  linear  feedback  of  the  state  variables.  The 
controls  are  modelled  as  linear  combinations  of  the  state  variables, 
but  the  machine  equations  are  not  linearized.  Therefore,  the  model 
parameters  which  are  assumed  to  be  time-independent  will  generally 
depend  on  the  strength  and  duration  of  the  pulse  contrary  to  the 
linearized  machine  models  where  the  equations  themselves  are  linearized. 


' 


Power  System 


Problem 


Steady-State  \ 


Dynamic 


— y 


Normal 

|  ; 

Faulted  I 

1  !  ! 

Normal 

j 

_ *  * 

Faulted 


Balanced 


Unbalanced 


1.  Algebraic  equations 

2.  System  assumed  linear 

3.  Instantaneous  photograph 
in  time. 


1.  Differential  equations 

2.  System  nonlinear 

3.  Solution  in  time  domain 


Fig.  1.1  An  Organization  of  Power  System  Problems 
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This  linear  feedback  model  is  considered  because  it  leads  to  a 
simple  implementation  compared  to  the  open  loop  system.  Due  to 
the  restrictions  on  the  form  of  the  control  as  mentioned  above, 
additional  constraints  have  to  be  considered  in  the  implementations 
of  the  optimal  solution  as  compared  with  the  unrealistic  model. 

The  method  of  solution  is  to  cast  the  system  into  the 
Pontryagin ' s  minimum  principle.  This  ensures  that  the  necessary 
and  sufficient  conditions  for  optimality  are  met  [1]  which  is 
important  for  obtaining  the  numerical  solution  for  such  a  highly 
nonlinear  system.  The  nonlinear  system  equations  and  the  constrains 
on  the  feedback  are  augmented  to  the  cost  functional  which  is  in  the 
quadratic  form. 

Then  using  Pontryagin 1 s  minimum  principle  and  the  gradient  descent 
method,  the  optimal  feedback  parameters  are  obtained  by  solving  a 
two  point  nonlinear  boundary  value  problem  (T.P.N.B.V.P. ) .  The  values  obtained 
for  these  parameters  depend  on  the  strength  and  duration  of  the 
disturbance  since  the  model  is  nonlinear  contrary  to  the  usual 
feedback  control  of  a  linear  model  as  shown  in  Fig.  1.1 

Using  this  kind  of  feedback,  ill-conditioning  is  handled  more 
easily  when  the  numerical  value  of  the  optimal  solution  is 


obtained . 
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CHAPTER  II 


BACKGROUND  OF  OPTIMAL  CONTROL 
2 .1  Introduction 

In  the  classical  control  system  design, acceptable  performance 
is  generally  defined  in  terms  of  time  domain  criteria  such  as  rise 
time,  settling  time,  and  peak  over shoot ; and  frequency  domain  criteria 
such  as  gain  margin,  phase  margin  and  bandwidth.  The  tools  used  are 
frequency  response,  root  locus  phase  plane  and  describing  function, 
although  these  methods  are  widely  used,  they  can  be  applied  only  to 
t ime- invar iant  systems,  with  unconstrained  variables. 

Modern  systems  are  multiple  input-multiple  output;  constraints 
on  state  and  control  variables  have  to  be  imposed.  The  objective 
of  optimal  control  theory  can  be  stated  as  follows:  Determine  the 
control  signals  which  will  enable  a  process  or  a  plant  to  satisfy 
certain  specified  physical  constraints  and  at  the  same  time  optimize 
some  index  or  measure  of  performance. 


controls  u(t)„ 


— > 


output  _y_(t) 


states  x(t) 


Figure  2.1.1  Diagram  of  Control  System 


5 


6 


The  optimization  problem  can  be  expressed  as: 

Given  a  dynamical  system  as  shown  in  Fig.  2.1.1, 

State  equation:  2i(t)  =  f_[x(t)  ,  u(t),t] 

Output  equat  ion  :^_(t)  =  g[x(t),  u_(t),t] 

with  (a)  the  initial  state  expressed  as  x(t  )  at  t=t 

~  o  o 

(b)  the  final  state  expressed  as  S(x(t^))  where  S  is 
called  the  target  S 

(c)  constraints  on  control  variables  expressed  as  u_(t)  e  U 

(d)  a  performance  criterion  or  cost  functional  usually  written  as 
J=J[x(tQ),  tQ,  u] 

=  /  L[x(t) ,  u(t) , t]dt 

t 

o 

2 . 2  Calculus  of  Variations 

2.2.1  Basic  Notions  of  Calculus  of  Variations 
PROBLEM 

Determine  the  extremal  x*(t)  for  the  functional 


J  =  /  L(x( t) ,  x(t) ,  t) dt 

t 

o 


(2. 2. 1.1) 


where  x(t)  is  a  vector 

t  and  t£  are  specified, 
o  t 

(1)  Let  x(t  )  =  x  and  x(t  )  =  x  . 

o  o  ft 

(2)  Assume  that  x(t)  and  x(t)  are  continuous. 

(3)  The  integrand  L  has  continuous  first  and  second  partial 
derivatives  with  respect  to  its  arguments. 

Consider  an  arbitrary  function  n(t) .  Assume  that  q  and  n(t) 
are  continuous  in  the  closed  interval  [tD>  ^  ]. 
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Let 

x(t)  =  x*(t)  +  e  n(t)  (2. 2. 1.2) 

where  e  ^0  is  a  real  parameter  and  x*(t)  is  the  solution. 

Substitute  (2. 2. 1.2)  into  (2. 2. 1.1) 

J(x)  =  J[x*(t)  +  e  n(t)] 

=  J(e) 
t 

=  /  L[t,  x*(t)  +  e  n(t) ,  x*(t)  +  e  n(t)]dt 

t 

o 

Since  x*(t)  extremizes  J,  we  can  say  that  J  is  a  min.  or  max. 
when  x(t)  =  x*(t)  or  when  e  =0  .  Considering  J  as  a  function  of  e 
we  set 


dJ 

de 


=  0 


£=0 


to  obtain  the  necessary  conditions  for  x*(t)  to  be  an  extremal.  That  is 


=  /  f [L  ^  +  L  . ,  ‘  n  +  L.  .  n]dt  =  0 

de  tJ  t  de  x*+eri  x>c+eri 


Since  =  0, t  does  not  depend  on  e.  The  above  can  be  simplified 
to  yield  at  e=0, 


dJ 


.  =  /  f  L  ,  •  n  dt  +  J  L.  •  n  dt  =  0 

de  tJ  x*  x* 


(2. 2. 1.3) 


If  the  second  integral  on  the  right  hand  side  is  integrated  by  parts 


'■V 
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and  substituted  into  (2. 2. 1.3),  We  obtain 


dJ  = 
de  t 


t 


t 


o 


f 


0  (2. 2. 1.4) 


Since  n(t)  is  an  arbitrary  function,  equation  (2. 2. 1.4)  will  be 
satisfied  only  if  the  following  necessary  conditions  are 
satisfied : 


o 


These  necessary  conditions  are  shown  in  Fig.  2.2.1. 

Since  x(t)  =  x*(t)  +  en(t)  for  all  t  in  [t  ,  t£]  we  have 

o  f 

x(t  )  =  x*(t  )  +  E  T1  (t  )  and 
o  o  o 

x(tf)  =  x*(tf)  +  e  T1  (t  )  . 

2.2.2  Variational  Notation 

We  replace  x(t)  by  x(t)  +  en(t)  and  represent  en(t)  by  6x(t). 
a.  First  Variation  of  a  Functional  L(x,x,t) 

By  definition 

6L  =  L(x  +  en,  x  +  ef),  t)  -  L(x,  i,  t) 

Expanding  L(x  +  ep,  x  +  ep,  t)  in  a  McLaurin's  series  in  e  and  retaining 
only  the  linear  terms; 


. 


. 


x(t)ik 


Fig.  2.2.1  Curves  of  Fixed  and  Free  End  Points 
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L(x  +  en,  x  +  en,  t)  =  L( 


x,  x, 


t)  + 


9L 


9(x+en)  de 


—  (x  +  en) 


e=o 


+ 


3L _ 

9(x+en) 


3  (x+en) 

9e 


e=0 


•  £ 


+  higher  order  items. 

Therefore,  6L  =  L  *(6x)  +  L.  (6x) 

x  x 

b.  First  Variation  of  J 


6J  =  J(x  +  en,  x  +  En,  t)  -  J(x,  x,  t) 


=  J  f  [L(x  +  En,  X  +  En,  t)  -  L(x,  x,  t) ] dt 


=  J  £  (SL)dt 


=  J  £  [L  (6x)+L.  (6x)  ]dt 

t  X  X 


=  /  f  L  (6x)dt  +  /  f  L.(6x)dt 

tJ  x  tJ  x 


Zf  ft,*..* 


=  tJ  Lx(6x)dt  +L^(6x)|t  -J  (6x)(—  L^)dt 


o  o 


A 

=  f  6x[L  -  L.  ]dt  +  L.(6x) 
tJ  x  dt  x  x  '  t 

o  ° 


\!  f  6xtLx  -  5F  Lx]dt  +  Lx(tf)6x(tf)  -  Li(to)4x(to) 


The  necessary conditions 


for  x(t)  to  be  extremal  (6J=0)  of  J  are 


- 


=  0 


2 .3  Pontryagin's  Minimum  Principle 

Pontryag in  ’  s  minimum  principle  is  an  extension  of  the  methods 
and  results  of  variational  calculus  to  the  case  of  bounded  control 
variables.  If  the  set  of  admissible  controls  is  unrestricted,  the 
calculus  of  variations  can  be  used  to  derive  necessary  conditions 
which  characterize  the  optimal  solution.  When  the  admissible  control 
set  is  bounded  in  some  way,  unrestricted  variations  in  u(t)  are  not 
allowed . 

The  minimum  principle  can  be  considered  as  an  extension  or 
generalization  of  the  calculus  of  variations  to  enable  one  to  take 
account  of  systems  whose  input  signals  have  constraints  of  certain 
types.  Consider  a  system  described  by  the  state  equation 


x  =  f (x ,  u,  t) 


(2.3.1) 


with 


x  ( t  )  =  x 
o 


o 


and 


x(tf)  =  xf 


x  is  the  n-dim.  state  vector. 


u  is  the  m-dim.  control  vector. 


The  cost  function  to  be  extremized  is 


■ 


NS 
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J  =  G[x(tf),  t  ]  +  /  f  L(x,  u,  t)dt  (2.3.2) 

t 

o 

J,  G  and  L  are  scalar  functions  of  their  arguments. 

Consider  the  Lagrange  problem  in  which  G=0 .  Introduce  a  vector  X ( t ) 
which  is  the  costate  vector  and  rewrite  (2.3.2)  as  follows: 

T 

J  =  J  [L  +  A  (f  -  x)]dt  (2.3.3) 

t 

o 

By  replacing  (2.3.2)  by  (2.3.3),  we  have  converted  a  constrained 
minimization  problem  into  an  unconstrained  minimization  problem. 

The  variation  6J  in  J  is: 


t 

o 


Since  6x(t^) ,  6x  and  6u  are  arbitrary  variations  the  requirement 
that  6J  becomes  zero  yields  the  following  conditions  for  x,  u_  to  be 
optimal  for  x,  u_  to  minimize  J. 

Setting  the  coefficient  of  6x  equal  to  zero  we  obtain 


<S>T  *  >’  €>  * » 


or 


(^T  A 
■(8x)  A 


3L 

8x 


(2. 3. A) 


Setting  the  coefficient  of  6u  equal  to  zero  we  obtain 


(it )T 


=  0 


(2.3.5) 


•  '  >*  (7T>  ‘  -  ^'))y  1  td 
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A  Hamiltonian  function  is  defined  as  follows 


H(x,  A,  u,  t)  =  L(x,  u,  t)  +  AT  f(x,  u,  t) 


Differentiating  (2.3.6)  with  respect  to  A 


®  -  f 
3  A 


=  x 


or 


x  = 


_3H 

3  A 


Differentiating  (2.3.6)  with  respect  to  x 


_3H  _  3L  ,3f.  T  . 

3x  3x  3x 


Comparing  (2.3.4)  and  (2.3.7),  we  obtain 


A  =  - 


_3H 

3x 


Differentiating  (2.3.6)  with  respect  to  u 


3H  _  3L  ,3f.T  . 

^ +  W  A 


Comparing  (2.3.5)  and  (2.3.8).  We  obtain 


(2.3.6) 


(2.3.7) 


(2.3.8) 


.  .  ! 
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2.3.1  Continuous  Optimal  Control  Problem 
Free  end  point  case 


x(t)  A 

I 

I 

I  !  I 

I  ! 

x(t^)  is  free 

. 


o 


t 

o 


Fig.  2.3.1  Free  End  Point  Case 

Let  u(t)  be  an  admissible  control  and  x(t)  be  the 
corresponding  trajectory  of  the  system  described  by 

x(t)  =  f(x,  u,  t) 

Let  x(tQ)  =  xq,  tQ  and  tf  be  specified  while  x(tf)  is  free  as 
shown  in  Fig.  2.3.1. 

The  necessary  condition  (a)  and  sufficient  condition  (b)  for 

u(t)  to  be  the  optimal  control  which  takes  x(t)  from 

x  to  some  state  xr  at  tr  while  minimizing  the  function  J 
o  — f  f 


15 


t 

J  =  I 


f 


t 

o 


rp 

[x(t)  Q  x ( t )  +  u(t)T 


R  u (t) ] dt 


are: 


(a)  There  exists  a  function  or  vector  A(t)  such  that  x(t)  and  A(t) 
are  solutions  of  equations 


x 

and 


m 

3  A 


A  = 


3H 

3x 


and 


subject  to  the  boundary  conditions 


x(t  )  =  x 
o  o 

A(tf )  =  Af 

(b)  The  matrix  Q  must  be  positive  definite  or  positive  semidefinite 
for  t  in  [tQ,  t^]  in  order  to  establish  the  sufficient  condition 
for  a  minimum.  The  matrix  R  must  also  be  positive  definite 
for  t  in  [ t  ,  t  ] . 

O  1  A 

The  functional  H(x,  A,  u,  t)  has  a  minimum  if  — tt  >  0  and  —  =  0. 

3u2  3u 

2 .4  Direct  Computational  Method 
2.4.1  Gradient  Techniques 


Optimization  methods  can  be  divided  into  direct  and  indirect 
methods.  The  first  one  comprises  various  gradient  methods  whereby 


. 
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the  gradient  of  the  cost  functional  with  respect  to  the  control 
function  is  used  to  successively  estimate  the  optimal  control.  The 
second  one  is  the  quasilinearization  method  which  involves  two  point 
boundary  value  problem  by  a  sequence  of  solutions  of  a  system  of  linear 
differential  equations.  This  is  denoted  as  a  direct  method  of 
computation  of  optimal  control. 

Consider  a  point  x  for  which  H(x)  is  not  a  relative  minimum. 

The  Taylor  series  expansion  is  known: 

H(x  +  h)  =  H(x)  +  hH'(x)  +  R? 

where  the  remainder  terms  are  . 

Choose  h  =  -  aH’(x),  a>o,  so  that: 

H(x  +  h)  =  H(x)  -  a[H’(x)]^  + 

Now,  if  a  is  sufficiently  small  such  that 

|aH’(x)|<e,  then 
| R2 |  <  a[H* (x) ]2 

and 

H(x  +  h)  <  H(x) . 


This  clearly  indicates  therefore  that  if  we  start  with  an 
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arbitrary  initial  guess  x°  for  x*,  where  H(x*)  is  a  relative 
minimum  and  update  each  successive  estimate  x1  according  to  the 
iterative  algorithm: 


x 


(2. 4. 1.1) 


then 


H(x1+1)  <  HCx1) 


and  x1  will  eventually  converge  to  x*.  The  vector  H’(x)  represents 


the  gradient  or  the  direction  of  steepest  descent.  Equation  (2. 4. 1.1) 
is  the  gradient  algorithm. 

2.4.2  Direct  Computation  of  Optimal  Control 

We  assume  that  there  are  no  inequality  constraints,  that  the 
initial  and  terminal  times  are  fixed,  and  that  the  initial  state  is 
fixed  and  the  terminal  state  is  unspecified. 

Thus  we  wish  to  minimize 


t 


t]dt 


(2. 4. 2.1) 


o 


for  the  system 


x  =  f (x,  u,  t) ,  x(tQ)  =  xc 


(2. 4. 2. 2) 


We  define  the  Hamiltonian 


T 

H(x,  A,  u,  t)  =  L(x,  u,  t)  +  A  f(x,  u,  t) 


(2. 4. 2. 3) 


8  * 


■ 
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and  then  set 


3H 

3x 


x  =  -3-LJ-^u't)  +  MO  [-3-fT(3-u>t)] 


and  then  set 


-  A  =  +  A(t)  [itil^O.] 


(2. 4. 2. 4) 


with  the  terminal  condition 


A(tf) 


3G[x(tf) ,tf ] 
3x(tf) 


(2. 4. 2. 5) 


Also,  we  minimize  the  Hamiltonian  with  respect  to  a  choice  of  u. 
For  the  case  where  any  control  is  admissible,  we  use 


3H 

3u 


0  = 


3L(x ,u , t) 
3u 


+  A(t) 


T 

3f  (x,u,t) 
3u 


(2. 4. 2. 6) 


We  shall  now  guess  a  solution  for  the  optimal  control,  and 

3H 

therefore,  we  will  not  obtain  —  =0.  We  will  solve  the  differential 

3u 

system  equality  constraint  (2. 4. 2. 2)  forwards  in  time  from  t^  to 
t^  with  the  assumed  value  of  u  and  will  also  solve  the  adjoint 
equations  (2. 4. 2. 4)  backwards  in  time  from  t^  to  t  with  the  terminal 
conditions  (2. 4. 2. 5). 

Thus,  the  incremental  first-order  change  in  the  cost  function 
(2. 4. 2.1)  becomes  for  a  control  differing  by  an  amount  6u(t) 


from  u(t) 


' 
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«J  -  /  f  «u(t)dt 


(2. 4. 2. 7) 


If  we  wish  to  make  the  largest  change  in  6j,  we  would  calculate 
3H 

the  gradient  —  and  then  make  6u  directed  opposite  to  the  gradient 


Su(t)  =  -  a 


9u(t) 


(2. 4. 2. 8) 


The  change  is  such  as  to  make  J  smaller  if  a  >  0,  which  is  desirable 
since  we  wish  to  minimize  J.  To  start  the  procedure,  we  assume 
that  we  have  some  nonoptimal  control  u1(t) ,  we  determine  x1(t)  by 
solving  (2. 4. 2. 2)  and  A(t)  by  solving  (2. 4. 2. 4)  with  the  terminal 
condition . 

The  predicted  change  in  J,  6J1,  may  if  desired  be  calculated 
from  (2. 4. 2. 7).  A  new  trial  value  of  u 


ui+1(t)  =  ui(t)  +  6ui(t) 


is  then  found  and  the  procedure  repeated  until  either  the  control  or 
the  cost  function  does  not  change  significantly  from  iteration  to 


iteration. 


■ 


OJ 


■ 


CHAPTER  III 


MODELLING  OF  THE  SYNCHRONOUS  MACHINE 
3.1  Park's  Transformation 


q  axis 


a 


d  axis 


a  axis 


Fig.  3.1  Pictorial  Representation  of  a  Synchronous  Machine 

As  shown  in  Fig.  3.1,  the  d  axis  of  the  rotor  is  defined  at 
some  instant  of  time  to  be  at  angle  0  rad  with  respect  to  a  fixed 
reference  position.  Let  the  stator  phase  current  i^,  i^  and  i  be 
the  currents  entering  the  generator  terminals.  If  these  currents 
are  projected  along  the  d  and  q  axis  of  the  rotor,  the  relations 
are  obtained  as: 


20 


21 


2 

i,  =  t  [i  cos  e  +  i  cos(e  -  120°)  +  i  cos'e  +  120°)] 
d  3  a  b  c 

2 

i  =  t  [i  sin  0  +  i  sin(e  -  120°)  +  i  sin(e  +  120°)] 
q  o  a  b  c 


Note  that  for  convenience  the  axis  of  phase  "a"  was  chosen  to  be 
the  reference  position  [  15 ] . 

The  effect  of  Park's  transformation  is  simply  to  transform  all 
stator  quantities  from  a,  b,  and  c  into  new  variables  the  frame  of 
reference  of  which  moves  with  the  rotor.  Park’s  transformation 
uses  two  of  the  new  variables  as  the  d  and  q  axes  components.  The 
third  variable  is  a  stationary  current,  which  is  proportional  to 
the  zero-sequence  current. 

This  is  given  by  [  15] 


l  = 


1_ 

3 


(i  + 


+  V 


where  i  produces  no  space  fundamental  field  in  the  air  gap.  Thus, 
o 

by  definition. 


* 


\ 


i 

o 


i 

q 


where  the  Park's  transformation  P  is  defined  as  [l8]: 


■ 


4 
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N 


P 


cos  0  cos( 0-120°)  cos( 0+120°) 


(3.1.1) 


I 

sin  0  sin( 0-120°)  sin( 0+120°) 

✓ 


If  the  transformation  (3.1.1)  is  unique,  an  inverse  transformation 
also  exists  wherein  one  may  write 


f  J  "1 

1 

•  'i 

1 

a 

O 

i. 

i— i 

i 

Pm 

(1 

i , 

b 

d 

i 

. 

i 

c 

<  q  J 

where  the  inverse  may  be  computed  to  be 

c-os  0  sin  0 

cos( 0-120°)  sin( 0-120°) 

cos(e+120°)  sin ( 0+120°) j 

Consider  the  expressions  for  the  flux  linking  the  coils  shown 
in  Fig.  3.1.  These  are  given  in  matrix  form  in  equation  (3.1.2). 


In  this  equation. 


. 

■ 
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/ - 

cd 

r< 

L 

aa 

Lab 

L 

ac 

1 

1 

LaF 

LaD 

LaQ 

1 

f  i 

a 

stator 

\ 

= 

L, 

ba 

Lbb 

Lbc 

t 

1 

I 

LbF 

LbD 

O' 

*b 

X 

c 

L 

ca 

Lcb 

L 

cc 

1 

i 

i 

i 

LcF 

LcD 

LcQ 

i 

c 

i 

X 

F 

LFa 

hb 

LFc 

i 

1 

l 

1 

L 

FF 

lfd 

lfq 

i 

rotor 

X 

D 

LDa 

LDb 

LDc 

i 

i 

i 

ldf 

ldd 

ldq 

h 

\ 

O' 

r< 

lLQa 

LQb 

LQc 

1 

1 

i 

1 

lqf 

lqd 

LQqj 

/q  ; 

stator 


rotor 


3.1.1  Self  Inductance  of  Stator 

The  phase-winding  self-inductances  are  given  by: 


(3.1.2) 


=  A 

+  A0 

cos2q 

aa  o 

2 

i. .  A 

+  A„ 

cos2 ( 0-120°) 

bb  o 

2 

=  A 

cc  o 

+  a2 

cos2(0+12O°) 

where  A  and  are  constants, 
o  2 

3.1.2  Mutual  Inductances  of  Stator 

The  phase -to-phase  mutual  inductance  are  functions  of  Q  but 
are  symmetric. 


ab 

Lba  - 

-  B 

o 

B2 

cos2 ( 6+30°) 

'be 

Lcb  = 

-  B 

o 

B2 

cos2 ( e-90°) 

ca 

Lac  = 

-  B 

o 

-B2 

cos2 ( 0+150°) 

where  B  and  Bn  are  constants, 
o  2 
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3.1.3  Self-Inductances  of  Rotor 

Since  saturation  and  slot  effect  are  neglected,  all  rotor  self¬ 
inductances  are  constants. 


3.1.4  Mutual  Inductances  of  Rotor 

The  mutual  inductance  between  windings  F  and  D  is  constant  and 
does  not  vary  with  9.  The  coefficient  of  coupling  between  the  d  and 
q  axes  is  zero  and  all  pairs  of  windings  with  90°  displacement  have 
zero  mutual  inductance.  Thus, 


0 

0 


3.1.5  Mutual  Inductances  Between  Stator  and  Rotor 

All  of  these  are  functions  of  the  rotor  angle  9.  From  the 
phase  windings  to  the  field  winding. 


L  =  L  ,  =  C  cos(9-120°) 
bF  Fb  1 

L  =  U,  =  C-  cos(9+120°) 
cF  Fc  1 


' 
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In  what  follows  the  damper  windings  are  neglected,  however, 
the  effect  due  to  these  windings  is  taken  into  account  by 
introducing  a  damping  factor  in  the  dynamical  equation.  See 
equation  (3.3.1). 


3 .2  Flux  Linkage  Equation 

The  sinusoidal  flux  wave  exists  in  the  air  gap  and  it  has  a 
magnitude  and  angle  in  terms  of  this.  One  obtains  more  insight  into  the 
constants  of  the  machine  when  one  resolves  this  flux  into  the 
components  q,  d  on  the  rotor  and  one  can  write  the  voltage  E  in  terms 
of  these  components  as  follows  [26]: 

E  =  p ( .cos0  +  \p  sine)  +  (R  +  L-p)i  -Lpi,  -  Lpi 
a  r  md  mq  a  l  a  m  b  m  e 

where  p  is  ^7 

is  the  leakage  inductance  of  coil  "a" 

L  is  the  part  of  mutual  inductance  between  two  armature  coils 
m 

due  to  flux  that  does  not  cross  air  gap. 

Using  i  4-  i,  +  i  =  3i  and  i  in  terms  of  i  ,  i  ,  one  obtains, 
°abc°  a  dq 


'  • 


. 
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Ea  =  p[(Vd  +  Laid)  COSe  +  +  L-i-)  sin0] 


mq  a  q 


+  L  pi  +  R  i 
o  o  a  a 


=  p  [\J>  cose  +  ^  sine]  +  L  i  +  R  i 
a  q  o  o  a  a 


where 


L  = 


L_  +  L 
1  m 


L  = 

L  • 

-  3L 

o 

a 

m 

^d  = 

*»& 

+  L  i, 
a  d 

i|j  = 

+  L  i 

q 

mq 

a  q 

\p , ,  ip  are  total  flux  linkages  with  an  armature  coil  located  on  the 
d  q 

appropriate  axis  due  to  both  the  main  air  gap  flux  and  the  armature 

leakage  flux.  is  the  effective  leakage  inductance  of  either  of 

the  axis  coils.  L  is  the  inductance  associated  with  the  zero- 

o 

sequence  current . 

If  one  substitutes  for  E  ,  i  in  terms  of  E , ,  E  ,  E  ,  i , ,  i  ,  i  , 

a’  a  d  q  o  d  q  o 

one  gets 


(Ed  -  Raid)  cos  0  +  (E^  -  Raiq)  sin  6  +  (EQ  -  Raip  -  LQp  iQ) 

=  cos  0(piK  +0)Tp  )  +  sin  6(p^  -u)  iK) 

da  q  a 


where  p  9  =  w. 


.  '  I 
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Equating  coefficients  of  sing,  cos0  terms  independently,  one  gets 


Ed  =  p*d  +  u.*q  +  Raid 

E  =  -  0)^,  +  pij;  +  R  i 
q  d  q  a  q 

EL  =  (R  +  L  p)i 
a  o^  o 


(3.2.1) 


Now  one  must  relate  the  fluxes  to  the  current  i , ,  i  ,  i-  and 

d  q  d  q  f 

also  define  the  flux  \p  and  also  find  the  expression  for  V^. 


*d  =  Lfd  *f  +  Ld  Ad 


= 


L  i 

q  q 


V  ^  =  (R^  +  pLff)if  +  L^p  i 


fd1 


where 

L-,  is  mutual  inductance  between  coil  f  and  d 
f  d 

L,  is  self  inductance  of  coil  d 
d 

L  is  self  inductance  of  coil  q 

q 

The  above  equations  are  easier  to  handle  if  one  uses  per  unit  (p.u.). 
In  this  case 


Jf  d 


=  L 


md 


L  ,  +  L 
md  a 


Jff 


Lmd  +  Lf 


L  +  L 
mq  a 


where 


L  is  per  unit  mutual  inductance  on  the  d  axis 
md 

L  is  per  unit  mutual  inductance  on  the  q  axis 
mq 

L  ,L  are  leakage  per  unit  inductances, 
a  f 


* 
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^  =  (L  +  L  )i,  +  L  ,  i. 
d  md  a  d  md  f 

^  =  (L  +  L  )  i 

q  mq  a  q 

Vf  =  LmdP  *d  +  [Rf  +(Lmd  +  Lf)pUf 


(3.2.2) 


Let  ♦f  "  Lmd  +  (Lmd  +  LfHf>  then  vf  =  Rf  *f  +  P  V  FrOT1 
equation  (3.2.2),  the  flux  linkage  equation  is  obtained  as 


d_ 

dt 


(uj  ljO  =  to  V. 
of  of 


CO  R  i 
off 


(3.2.3) 


3 . 3  Torque  Equation 


M  + 
e 


doo 

J  dF  +  Kd 


0) 


(3.3.1) 


where 

is  the  total  instantaneous  applied  torque 

is  the  total  instantaneous  electrical  torque 

K,co  is  the  damping  torque 
d 

J  is  the  moment  of  inertia. 

is  positive  if  power  is  passing  into  the  machine  from  outside  at 

a  positive  speed.  If  the  machine  accelerates  constantly  from  the  rest 

to  co  in  one  second,  then 
o 


CO 

o 


where  the  inertia  constant  H  is: 


! 


.  • 
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stored  energy  at  synchronous  speed  in  KW-sec 

Rated  KVA 


Substituting 


into 


equation  (3.3.1), 


M 


t 


2H 

+  —  p  co  +  K , 

(jo  d 

o 


0). 


Multiplying  (3.2.2)  by  w  and  redefining  the  constants  in  terms 

o 

of  impedances  instead  of  inductances,  one  obtains: 


co  ip  =  x  i  +  x  i 
o  d  d  d  md  f 

w  \p  =  x  i 
o  q  q  q 

co  \p  =  x  i  +  x  i 
of  md  d  fd  f 


where 


x  ,  = 

0) 

(L 

d 

o 

md 

x  = 

0) 

(L 

q 

o 

mq 

Xmd= 

o 

3 

^md 

Xfd= 

0) 

O 

(Lmd 

x  „  +  x 
md  a 

x  +  x 
mq  a 


Xmd  +  Xf 


It  can  be  shown  that  if  the  machine  is  connected  to  an  infinite 


bus  by  a  series  impedance  x  ,  R  then  one  can  change 


x  =  x '  +  x 
a  a  e 

R  =  R'  +  R 
a  a  e 
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(The  primed  quantities  are  original  machine  parameters).  Consider 

the  machine  connected  directly  to  the  infinite  bus,  where  E,  and  E 

*  d  q 


Fig.  3.3  Single-Line  Diagram 

are  the  components  of  the  infinite  bus  voltage  in  the  appropriate  axis. 
Now  it  is  assumed  in  the  voltage  equations  that  as  shown  in  Fig.  3.3 

Pl>d  =0 

p%  =0 

R  =0 
a 

R  =0 
e 

03  =  0)  . 

o 

Thus,  it  follows  that 


0)  \p 
o 

03 

O 

03  \b 

o 


d 

q 

f 


Xd  h  +  Xmd 


X  1 

q  q 


Xmd  h  +  Xfd 


(3.3.2) 


■ 
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where 


X,  =  x'. 

+ 

X 

d  d 

e 

X 

II 

X 

+ 

X 

q  q 

e 

(The  primed  quantities  are  original  machine 
parameters) . 


When  the  bus  voltage  is  substituted  in  these  equations,  two 
dynamical  equations  are  obtained  as  shown  below.  Let  the  voltage 
in  phase  "a"  be 


E  =  E  Pi  sin  a)  t  =  E  sin  oo  t 
a  rms  o  m  o 

where 

sin  a)  t  =  sin(0-b  <5)  =  sin  0  cos  6 cos  6  sin  6 
o 

then 

E  =  E  sin  0  cos  6  4-  E  cos  0  sin  6  (3.3.3) 

am  m 

Now  from  inverse  Park's  transformation,  one  has:  [Z(-j 


E  =  cos  0  E ,  —  sin  0  E  (3.3.4) 

a  d  q 


Comparing  (3.3.3)  and  (3.3.4),  one  has 


E,  =  E  sin  6  =  Pi  E  sin  6 
d  m  rms 

E  =-E  cos  6  =  ~Pi  E  cos  6 
q  m  rms 


Now  one  expresses  the  currents  id,  i  ,  if  and  the  fluxes 


„ai3 
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u)  ip, ,  “  ^  in  terms  of  the  bus  E  ,  the  angle  6  and  w 
o  a  o  q  rms*  b  of 


a)  =  /2  E  cos  6 

o  d  ms 

w  ip  =  Jl  E  sin  6 
o  q  ms 


1  = 

q 


/2  E  sin  6 
ms 


x 


From  equation  (3.3.2), 


Xd  I 


-1 


Xd  Xmd 


L Xmd  Xf d  J 


Jl  E  cos  6 
ras 

“o*f 


X 


fd 


XdXfd  Xmd 


md 


x 


md 


s  xdXfd  Xmd 


X  dXfd  Xmd 


x 


XdXfd  Xmd 


ll  E  cos  6 
rms 


} 


“o*£ 


Therefore 


1 ,  = 


Jl  E 


rms 


x 


2  cos6- 


x 


md 


md 


x 


fd 


xj - 

d  x 


fd 


-)  u>  ipJ 


x ,  - 


x2  7  0Tf 
md 


x 


fd 


-x 


i„  = 


md 

"fd 


x 


x  - 

d  x 


— ~ — )  12  E  cos  6  H — 
x^,  ras  x 

md 


fd 


( - J5-)  »0*f 


fd 


x ,  - 
d  x 


md 


fd 


' 
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i  = 

q 


G  E  sin6 
rms 


x 


(3.3.5) 


co  Tp ,  =  'll  E  cos  6 
o  d  rms 

co  ip  =  v^2E  sin  6 
o  q  rms 


Substituting  (3.3.5)  into  M  =  \  (co  ifj.  i  -  oo  ip  i.)  ,  it  becomes 

0  2  od  q  o  q  d 


M  = 
e 


S.  sin  6  cos  6  +  Sr  co  \li_  sin  6 
4  5  of 


where 


(3.3.6) 


x 


S4  ’ 


E2  (x  -x  +  -^-) 
rms  q  d  x  ^ 


(  Xmd\ 
x  (x  -  - ) 

q  d  x. 


fd 


E  x  , 

„  _  rms  md  ,  1  v 

5  =  xfd  T7^ 

d  ^ 


Now  from  the  definition  of  6,  one  has  9  =  co  t  +  6 

o 

•  •  •  •  • 

0  =  6  =  co 

where  coq  is  the  synchronous  speed 

6  is  the  rotor  angle  with  respect  to  a  synchronous  reference, 


Mechanically  one  has 


M  ^-4-  +  K.  4^  =  M'  -  M 


dt 


2  d  dt  t  e 


where  6  is  the  rotor  position 


M  is  the  moment  of  inertia 
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K  is  the  damping 


is  the  mechanical  torque  of  prime  mover 


Mg  is  the  electrical  torque 


This  corresponds  to  the  electrical  phase  angle  of  the  internal 


voltage 


M 


d26 

dt2 


+  K 


d5 

dt 


=  M 


-  M 


e 


where  M  =  M’ 
t  t 


So  it  becomes 


K,  a)  is  the  net  prime  mover  torque, 
d  o 


(3.3.7) 


Substituting  (3.3.6)  into  (3.3.7),  it  becomes: 


=  S  [Mt  '  S4  SinS  C0S6  '  VVM  sin,S  '  Kd  ft] 

dt 


where 


a) 

o 


Substituting  (3.3.5)  into  (3.2.3),  it  becomes: 


d  (a)  ^  ) 

- =  a)  V,,  -  A (w  +  C  cos6 

dt  of  of 


where 
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x 


A  = 


x 


fd 


x 


)  0)  R 
o  f 


X  ,  - 

d  x 


md 


fd 


C  = 


/2  E  x  ,  03  Rr 
rms  md  of 


md. 

x  (x  -  - — ) 
f  d  d  x, 


fd 


CHAPTER  IV 


MATHEMATICAL  MODEL  OF  TURBO-GENERATOR 
4  *1  Machine  Equations  in  State  Space 

If  one  defines  6  =  rj ,  one  obtains  the  following  third  order  system 

from  the  previous  chapter. 


6  =  n 


U1  S4  S5  Kd 

h  =  Tf - 77-  sin  6  cos  5  -  —  x„  sin  6  -  — —  x_ 

MM  M3  M2 


x^  =  -  A  x3  +  C  cos  6  (4.1.1) 


Summarizing  the  assumptions  used  in  the  derivation  of  this 
third  system; 

1.  Damper  windings  were  not  considered,  instead  a  damping 
coefficient  is  introduced  in  the  torque  equation. 

2.  Transient  response  in  the  transmission  line  is  neglected. 

3.  The  resistances  in  the  stator  winding,  the  transformer  and  the 
transmission  line  are  neglected. 

4.  The  time  derivative  of  the  fluxes  ty.,  ip  is  neglected,  i.e. 

d  q 

•  • 

^ ,  =  0  ,  ^  =  0  . 

d  q 

5.  0)  is  approximated  to  in  the  voltage  equations. 

Equations  (4.1.1)  represent  a  synchronous  machine  to  a  very 

good  approximation  [11] .  When  use  is  made  of  the  transformations 


x^  =  sin  6 
x2  =  n 


(4.1.2) 


36 


37 


the  machine  equations  (4.1.1)  become 


X1  =  X2X4 


U1  S4  s5  Kd 

X2  =  M~  ~  >T  X1X4  "  X1X3  "  >T  X2 


-  A  x3  +  C 


where 


x.  =  1  -  x' 


(4.1.3) 


It  is  assumed  that  the  disturbance  is  modelled  by  introducing 

a  torque  pulse  in  the  right  hand  side  of  the  torque  equation  of 

magnitude  K  which  is  increased  in  steps  of  5%  of  steady  state 
s 

value  M^,  and  a  duration  t  (0.5  sec).  Hence  in  this  model,  the 
values  of  the  variables  at  t  =  t  are  the  steady  state  values. 

Therefore,  the  machine  equations  (4.1.3)  for  this  model  are 
modified  to 


X1  X2  X4 

*2  =  S  [ul  "  S4  xlx4  ‘  S5  X3X1  "  Kdx2  +  K[u(t'to)  '  u(t-T)] 

x„  =  u’  -  A  x0  +  C  x.  (4.1.4) 

3  2  3  4 


x, 


-  x2X;L 


. 
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where 


u 


2 


co  V 

o 


f 


Here  the  last  equation  in  (4.1.3)  which  is  algebraic  is  converted 
into  a  differential  equation  as  in  [2_]  in  order  that  the  equality 
constraints  are  unified  in  structure  so  that  it  is  easier  to  find  a 
unified  approach  for  the  implementation  of  a  numerical  solution  of 
the  system. 

4 . 2  Modelling  of  the  Turbine,  the  Governor  and  the  Exciter 

4.2.1  Steam  Turbine 

To  develop  the  dynamic  characteristic  of  a  steam  turbine, 
consideration  must  be  given  to  the  energies  associated  with  the 
steam  system  and  the  mechanical  rotating  system.  The  normal 
controlling  input  for  a  turbine  is  the  control  value  and  the  output 
is  the  mechanical  power  delivered  on  its  output  shaft. 

Usually  the  turbine  is  solidly  connected  to  its  generator  so  the 
inertial  effects  of  both  units  are  considered  together  as  additive. 
For  this  dynamic  model,  the  effects  of  boiler  pressure  changes  can 
be  ignored  so  that  the  steam  control  valve  can  be  considered  as  a 
power  controller. 

The  steam  turbines  in  an  electric  power  system  normally  consist 
of  two  or  more  steam  expansion  stages  with  moisture-separator  and 
reheat  stages  between  them.  The  transfer  function  relating  the 
power  developed  by  a  turbine  stage  to  the  steam  flow  of  the  stage 
is  given  by: [12] 


jt\  i  >  ,r  ■,  ,  ■  A  k  3  06  •5etrra  ™i  J»*afcUflOn 

■ 
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P 


P 


in 


00 


00 


s 


1 

1  +  T  (k)  s 

D 


where 

P^n(k)  is  the  power  (per  unit)  developed  by  the  k^  stage, 

P  (k)  is  the  flow  of  steam  (per  unit)  into  the  k^  stage, 
s 

t  ll 

T^(k)  is  the  characteristic  time  constant  of  the  k  stage. 


Fig.  4.2.1  Block  Diagram  of  Steam  Turbine 

For  this  steam  turbine  the  transfer  function  relating  input 
valve  movement  to  output  mechanical  power  is: 


—  =  - ± -  (4. 2. 1.1) 

s  1  +  T,s 

b 

Note  that  the  transfer  function  is  a  linearized  approximation 
based  on  conditions  at  a  specific  operating  level. 

The  development  of  the  differential  equations  for  a  typical 


double  expansion  turbine  is  shown  in  [12], 


The  same  techniques 
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can  be  employed  to  develop  transfer  function  (4. 2.1.1)  for  this  type 
of  unit. 


where 

is  the  mechanical  power 
Pg  is  the  steam  power 

is  the  time  constant  of  the  turbine 
when  use  is  made  of  the  transformations 


x 


5 


=  P. 
m 


x 


7 


=  P 


s 


the  turbine  equation  (4. 2.1.2)  becomes: 


(4. 2. 1.3) 


4.2.2  Speed  Governor 

The  speed  governor  used  on  the  generator  prime  mover  may  be 
considered  either  as  a  speed  controller  or  power  controller  depending 
on  the  mode  of  generator  operation.  If  the  electric  generator  is 
supplying  an  isolated  load,  the  governor  tends  to  operate  as  a  speed 
controller;  however  if  the  generator  is  connected  to  a  large  system 
where  the  total  generating  capacity  is  much  larger  than  that  of  the 
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+ 

->- 


1 

i 

Load 

Sensor 

Fig.  4. 2. 2.1  Functional  Elements  of  a  Speed  Governor 

unit  itself,  the  frequency  cannot  be  appreciable  altered  by  the 
individual  unit  and  therefore  the  governor  operates  as  a  power 
controller.  The  general  arrangement  of  the  functional  elements  of 
the  governor  is  shown  in  Fig.  4. 2. 2.1. 

For  speed  and  power  control  study  purposes,  both  the  steady- 
state  gain  and  the  transient  performance  of  governor  are  of  interest. 
To  develop  the  transfer  function  of  a  typical  speed  governor,  the 
steady-state  droop  characteristic  has  to  be  considered. 

The  steady-state  performance  of  a  governor  depends  on  the 
value  of  the  steady-state  speed  droop  R.  That  is,  under  steady-state 
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Fig.  4. 2. 2. 2  Governor  Droop 


condit ions 


where 

A V  is  the  amount  of  turbine  inlet  valve  movement 

R  is  the  governor  droop 

■jji-  is  referred  to  as  the  steady-state  gain 

Aw  is  the  speed  error. 

Fig.  4. 2. 2. 3  illustrates  the  most  common  transfer  function  for 
steam  turbine  governor  relating  speed  change  Aw  to  steam  control 
valve  movement  AV.  With  reference  to  the  governor  characteristic 
shown  in  Fig.  4. 2. 2. 2,  adjustment  of  R  in  Fig.  4. 2. 2. 3  has  the  effect 
of  adjusting  the  slope  of  the  characteristic  while  variation  of  w^ 
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AV 


Fig.  4. 2. 2.3  Block  Diagram  and  Transfer  Function  for 
a  Steam  Turbine  Governor 

adjusts  the  whole  characteristic  up  and  down  while  maintaining  the 
same  slope. 

From  Fig.  4. 2. 2. 3,  we  obtain: 

AV  ^  1/R  (4. 2. 2.1) 

A  co  1  +  T  s 
r  g 


where 

T  is  the  time  constant  of  the  governor 
g 

Td  is  A/R 
§■ 

A  is  the  governor  constant 

Substituting  P  ,  Y  into  (4. 2. 2.1)  V,  co  ,  the  transfer  function 
°  s  o  r 

(4. 2. 2.1)  becomes: 


G2G3 


1  +  T  s 
g 


(4. 2. 2. 2) 


where 
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p 

s 

is 

g3(v 

•  Mt) 

Y 

is 

to 

M 

o 

r 

t 

1/R 

is 

G2 

G2 

is 

the 

governor  gain  constant 

G3 

is 

the 

steam  valve  constant 

T 

g 

is 

the 

time  constant  of  the  governor 

Y 

n 

is 

the 

speeder  gear  setting. 

Y 


Fig.  4. 2. 2. 4  Block  Diagram  of  a  Steam  Turbine 
Governor 

The  transfer  function  (4. 2. 2. 2)  can  be  converted  into  the 
following  differential  equation. 


P 

s 


G2G3 


g 


1_ 

T 

g 


P 

s 


When  use  is  made  of  the  transformations 


(4. 2. 2. 3) 


P 

s 

Y 


o 
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the  governor  equation  (4. 2. 2. 3)  becomes 

*7  -  Y-  tG3G3ul  ■  x7]  (4. 2. 2. 4) 

g 

4.2.3  Exciter 

4 . 2 . 3 . 1  Automatic  Voltage  Regulator  (A . V . R . ) 

In  most  modern  generators  the  output  voltage  is  controlled  by 
automatic  devices  so  as  to  remain  at  a  constant  prearranged  value. 


Pilot  Exciter 


I 


Main  Exciter 

Fig.  4. 2. 3.1  Excitation  Arrangements  for  a 
Synchronous  Generator 


Rotor  Winding 


The  excitation  arrangements  for  the  rotor  of  a  synchronous 
generator  is  shown  in  Fig.  4. 2. 3.1.  There  are  two  broad  divisions 
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of  automatic  regulator  both  of  which  set  out  to  control  the  output 
voltage  of  the  synchronous  generator  by  controlling  the  exciter  field. 
In  general  the  deviation  of  the  terminal  voltage  from  a  prescribed 
value  is  passed  to  control  circuits  and  thus  the  field  current  is 
varied . 

One  type  can  be  broadly  classed  as  electromechanical.  Here  a 
voltage  proportional  to  the  deviation  voltage  operates  a  solenoid 
assembly  to  vary  the  pressure  exerted  on  a  carbon-pile  resistor  in 
the  exciter  field,  thus  varying  its  resistance. 


Fig.  4. 2. 3. 2  The  Effect  of  a  Dead  Band  in  Automatic 
Voltage  Regulator 

This  type  suffers  from  the  disadvantages  of  being  relatively  slow 
acting  and  processing  dead-bands,  i.e.  a  certain  deviation  must 

occur  before  the  mechanism  operates  as  shown  in  Fig.  4. 2. 3. 2. 

The  other  main  group  of  regulators  is  known  as  continuously 
acting  and  these  are  faster  than  the  above  and  have  no  dead  bands. 
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Deviation  Amplifier  Exciter  Generator 


1 

r  I . 

G2 

1 - 1 

_  c;  j  . 

G4 

G1  f 

i  3  1 

1  +  T2S 

Ll+.I3s  J— 

1  4-  T4  s 

1  +  T5s 
Feedback 

- * - - — * - 

Fig.  4. 2. 3. 3  Block  Diagram  of  a  Continuously 

Acting  Automatic  Voltage  Regulator 

A  general  block  diagram  of  a  typical  control  system  is  given  in 
Fig.  4. 2. 3. 3. 

4. 2. 3. 2  The  State  Equation  of  Exciter 

In  the  present  investigation,  an  exciter  is  assumed  to  be 
present,  but  there  is  no  feedback  applied  to  it  from  the  machine 
terminal . 

The  terminal  voltage  can  be  represented  only  as  a  complicated 
function  of  the  state  variables  6  and  and  cannot  easily 

be  included  in  the  dynamic  equations  of  the  system.  The  voltage  V 


ex 


. 
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applied  to  the  field  winding  of  the  exciter  is  assumed  to  be 
available  for  direct  manipulation  and  is  used  as  one  of  the  control 
input  variables. 

The  effect  of  the  exciter  is  approximately  represented  by 
a  single  time  constant  and  a  suitable  gain.  There  are  various  forms 
of  exciters  and  the  parameters  of  the  model  will  depend  on  the  type 
of  exciter  considered.  For  thyristor  rectifier  excitation  supplied 
from  transformers  or  a.c.  exciters  with  silicon  diode  rectified 
output,  the  time  constants  vary  between  0.05  and  0.3  sec.  while  for 
d.c.  exciters  larger  time  constants  of  1  to  2  sec.  are  common. 


Fig.  4. 2. 3. 4  Block  Diagram  of  Exciter 


The  transfer  function  of  the  exciter  is: 


C 

r 


1  +  T  s 
e 


03  V 
o  ex 


(4. 2. 3.1) 


9»J  ,n*V  3t>  lU 


- 
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where 


C  is  the  exciter  gain  constant 
r  b 


T  is  the  time  constant  of  the  exciter 
e 

to  V  is  the  exciter  field  voltage  times  the  synchronous  speed 
o  ex 

co £  is  the  field  voltage  times  the  synchronous  speed. 


In  the  time  domain  equation  (4. 2. 3.1)  becomes: 


co  V  =  —  [C  to  V  -  co  V  J 
of  T  roex  of 
e 


When  use  is  made  of  the  transformations 


x.  =  to  V 
6  of 

u_  =  co  V 
2  o  ex 


the  exciter  equation  (4. 2. 3. 2)  becomes 


x,  =  [C  u  -  xj 
6  T  r  2  6 


(4. 2. 3. 2) 


(4. 2. 3. 3) 


CHAPTER  V 


MATHEMATICAL  SOLUTION  OF  THE  MODEL 


5 . 1  System  Equations  in  State  Space 


Relays 

Servos 


Transmission 

Lines 

- W'. - — >-  Inf  in  it  e 

System 


Exciter 


Exciter  Field 
Voltage  V 

BX 

Fig.  5.1  Schematic  Diagram  of  the  System  Investigated 

The  model  used  here  is  a  specialized  form  of  the  mathematical 
model  (4.1.4)  derived  in  Chapter  3,  since  we  review  the  mathematical 
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model  in  Chapter  4  and  present  its  specialization  for  this  chapter. 

Inclusion  of  the  transfer  functions  for  the  turbine,  the 
governor  and  the  exciter  will  add  three  more  additional  differential 
equations  (4. 2. 1.3),  (4. 2. 2. 4),  (4. 2. 3. 3)  to  those  in  Section  4.2. 
Here  the  same  procedure  is  used  as  in  [12] .  The  system  equations 
for  a  machine  disturbed  by  a  torque  pulse  together  with  the 
differential  equation  generated  by  the  transfer  functions  are: 


X1  X2X4 


x. 


x. 


M 


M  X1X4 


S  K 

X  X  -  -p-  X  +  ^  [u(t-t  )-u(t-x)] 
M  3  1  M  2  M  o 


x„  =  x,  -  A  x_  +  C  x. 
3  6  3  4 


x4  =  -  X2X1 


(5.1.1) 


x5  =  tT  [X7  -  X5] 
b 


x6  =  t"  [Cru2  -  X6] 

e 


Xy  =  —  fG2G3ui  “  x7  ^ 


Here  x  ,  x2 ,  x^  and  x^  are  the  sin6,  derivative  of  torque  angle, 
field  flux  and  cos6  respectively,  x^,  x^  and  Xy  are  the  mechanical 
power,  the  field  voltage  times  the  synchronous  speed  and  the  steam 
power  respectively.  S^,  S,.,  K^,  M,  A,  and  C  are  the  machine  parameters 
derived  in  Chapter  3.  Tb»  T  ,  and  Te  are  the  time  constants  of  the 
turbine,  the  governor  and  the  exciter  respectively,  while  G^,  G^>  an<^ 


- 
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C^_  are  the  governor,  the  steam  valve,  and  the  exciter  gain  constants 
respectively. 

The  control  signal  is  assumed  to  control  directly  the  speeder 
gear  setting  of  the  governor  mechanism  and  thereby  the  valve  opening. 

One  of  the  control  input  variables  u^  is  assigned  to  Yq  which  represents 
the  speeder  gear  setting.  The  other  control  input  variable  u^  is 
the  voltage  V  applied  to  the  field  of  the  exciter  multiplied  by 

GX 

co  for  convenience  as  shown  in  Fig.  5.1. 
o 

Before  use  is  made  of  the  transformations 


sin6 

n 


co  ip 

o 

cos 


f 

6 


the  cost  function  to  be  minimized  is: 


.  f  s  2  2  ,  2  .  s, 

J  =  /  1  [a1(6-6  )  +  a2n  +  a3(x3"x3)  +  yi(x5~x5) 

t 

o 


+  y2(x6_x6)2  +  Y3(X7_X7^  +  Bl(urUir  +  62(u2'u2r] 


Sv  2 


sN  2 


sN  2. 


where  superscript  s  indicates  the  steady  state  value.  The  parameters 

ar  a  ,  a3,  y2,  Y3 ,  \  and  &2  are  penalty  factors  for 

deviations  from  steady  state  values. 


sN  2 


Due  to  the  transformation  (4.1.2)  one  should  express  (6  6  ) 
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g 

in  the  cost  function  in  terms  of  (x^-x^)  and  then  expand  the 
resulting  function  and  introduce  additional  pseudo-variables.  On 
the  other  hand,  one  could  alter  the  cost  function  slightly  and 
if  necessary  increase  the  factor  a,  to  compensate  for  this  change 
as  will  be  shown  below. 

s  2  s  2  s  2 

A  cost  function  with  (6-6  )  replaced  by  (x^-x^)  +  (x^-x^) 

is  mono tomic ally  increasing  in  the  range  (0,tt)  and  approximates 
s  2 

(6-6  )  to  a  high  degree  of  accuracy  over  the  range  (0,  tt/2)  .  This 
can  be  seen  as  follows: 


(x^x^)2  +  (x^-x^)2  =  2  [  1  —  cos(6-6S)] 


If  one  now  chooses  a  (a  =  2.5  a^)  for  the  coefficient  of 

[(x1-x®)2  +  (x^-x®)2]  in  the  cost  function  ,  one  will  have  the  same 

penalty  at  tt  for  both  the  old  and  new  cost  function  . 

It  should  be  noticed  that  in  this  case  small  errors  are  penalized 
more  heavily  than  before  which  is  a  beneficial  feature  of  the 
modified  cost  function.  Now  the  modified  cost  function  becomes. 


J  =  /  f  [o[(Xl-x®)2  + 
t 


,  s\2-\ 
(x4-x4)  ] 


+  a2x2  + 


(  S\2 
a3  X3  X3 


+  VX5_X5)2  +  T2(X6~X6^  +  Y3(x7"X7)  +  3l(ul"Ul) 


sx  2 


Sv  2 


Sx  2 


+  S2(u2-u®)2]dt 


(5.1.2) 


The  control  is  specialized  to  be,  for  this  problem,  of  the  form 
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s 

urui 


s 

u2-u2 


al(xrXl}  +  a2(x2_X2}  +  a3(x3"X3}  +  a4(x4-x4} 


^1^X1-X1^  +  ^2^X2~X2^  +  ^3^X3-X3^  +  ^4^x4_x4^ 


(5.1.3) 


where 


a^,  a^,  a^  and  a^  are  constants 
b^,  ^2*  ^3*  an<^  ^4  are  constants, 


However  the  magnitude  of  these  constants  in  general  depends  on 


the  strength  of  the  disturbance.  These  constants  clearly  satisfy 


5.2  Mathematical  Solution  of  the  Model 
By  substituting 


s 


yl 

xi 

X1 

CM 

= 

x2 

- 

s 

x2 

y3 

= 

X3 

- 

s 

X3 

y4 

= 

x4 

- 

s 

x4 

y5 

= 

X5 

- 

s 

X5 

y6 

= 

x6 

- 

s 

X6 

y7 

= 

x7 

- 

s 

X7 

V1 

= 

U1 

- 

s 

U1 

v2 

= 

u2 

— 

s 

u2 

(5.1.4) 


(5.2.1) 
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where 

y_^  and  are  deviations  from  the  steady  state  into  the  cost  function 
equation  (5.1.2)  and  equations  (5.1.1),  (5.1.3)  and  (5.1.4)  one 
obtains : 


t  f  f  r  r  2x  2  2  2  2  2 

J  =  J  [a(y;L+y4)  +  a2y2  +  a^3  +  y^  +  Y2y6  +  y3y? 

t 

(5.2.2) 

+  +  l  (v.a^  +  o.b^Jdt 

1=1 


Note  that  we  set  x2  =  0  because  the  torque  angle  does  not 

change  in  the  steady  state.  Also  quadratic  penalties  are  imposed 

on  the  a ! s  and  b ! s . 

l  l 


Yi  = 


X4y2  +  y2y4 


*2  =  I 


^5X3yl  ^<jy2 


S5Xly3 


S4Xiy4 


-  Vly4  '  S5yly3  +  K[u(t"to)  '  U(t_T)1 

y3  =  y6  '  A  y3  +  C  y4 

y4  =  '  yly2  ‘  y2Xl 

y5  =  T  Iy7  "  y5] 

b 

y6  ■  ?  [Cr  v2  -  y6] 
e 

y?  =  ~~  [G2G3  V;l  -  yy] 

g 


(5.2.3) 
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v.,  =  a. y.  4-  a_y_  +  a_y„  +  a.y. 
1  11  22  33  4'4 


v2  =  biyi  +  b2y2  +  b3y3  +  b 


a3  =  a_  =  a„  =  a.  =  0 
12  3  4 


bl  =  b2 


b  ,  =  b  =  0 
3  4 

Governor 


■ 

G2G3 

vi 

; 

T 

g 

— * — 

al  a2  a3  a4 


Turbine 


Generator 


yl  y2  y3  y4 


V2 


LOAD 


1  bi  b2  b3  bA 


Fig.  5.2  State-Space  Diagram  of  the  Model 
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It  should  be  noted  that  from  equations  (5.2.1)  and  (5.2.3)  it 
follows  that 

s  _ss  „ssn 
x5  -  S4X4X1  -  S5X1X3  =  ° 

xf  -  Axif  +  C  xf  =  0 

6  3  4 

at  steady  state. 

The  initial  values  are 

yl(to}  "  y2(to)  “  W  =  W  =  y5(t0)=  W  =  W  =  0 

(5.2.4) 


Then  the  augmented  cost  function  is  (Section  2.3): 


J1  =  /  f  +  a2y2  +  a3y3  +  yly5  +  Y2y6  +  V?  +  Vl 


+  B2v22  +  l  (v.a2  +  o.b2)  +  X]_(y1  -  [x°y2  +  y^]) 
i=l 


+  X2(y2  '  Sty5  "  S4x4yl  "  S5X3yl  "  Kdy2  '  S5xly3  "  S4xly4 


'  S4yly4  "  S5yly3  +  K(u(t-t0)  "  +  A3(y3_y6+Ay3_Cy4) 


g 

+  X^(y^  +  y  y  +  x1y2^  +  kl^Vl~^aiyl  +  a2y2  +  a3y3  +  a4y4 


)) 


+  +  b2y2  +  b3y3  +  b^))  +  ^(a^  +  ^(a^ 
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+  X?(a3)  +  Ag(a4)  +  AgCb^  +  A1Q(b2)  +  An(b3) 


+  X12(b4)  +  A13  [y  -  i-  (y7-y5)]  +  A14[y6  -  y-  (C  v  -y  >] 

b  e 


+  X15lh  -  T  (g2g3  vi-y7)]]dt 
g 


Integrating  by  parts  and  collecting  terms  which  are  quadratic 

and  linear  in  the  y’s  and  a!s  and  b ’ s  and  dropping  terms  which  do 

x  J  k 

not  depend  on  these,  the  augmented  cost  function  reduces  to: 


J  =  /  f[YTQY  +  V1RV  +  AXSA  +  BiTB  +  Ai(Y-F)]dt 


T  * 


where 


Y  =  [y1  y2  y3  y4  y5  y6  y?] 


v  = 


tvl  v2] 


A  =  ^al  a2  a3  a4 ^ 
bt  -  [b3  b2  b3  b4] 

T 

A  =  [A1  X2  X3  X4  X5  X6  A7  A8  A9  A10  A11  A12  A13  A14  A15] 


a 


a. 


0 


a. 


Q  = 


a 


Yi 


y* 


Y' 


0 
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R  = 


0 


N 


The  matrices  Q,R,S  and  T  are  required  to  be  positive  definite 
(Section  2.3)  for  minimum  . 

Thus,  the  Hamiltonian  is  defined  as  (Section  2.3): 


H  =  “(y^+y^)  +  “2y2  +  °3y3  +  yly5  +  y2y6  +  y3y7 


+  B1v^  +  +  Xi('[x4y2+y2y4]) 


+  V-  s  [y5  -  V**1  '  85X35,1  '  KdV2  '  85X15,3  '  84X35,4 


■ 


60 


"  S4yly4  "  S5yly3  +  K(u(t-tQ)  -  u ( t— t) ) ] ) 


+  A3(-y6  +  Ay3  -  Cy4)  +  \(y1y2  +  *[y2) 


+  k  (v,  -  (ay  +  a_y_  +  a_y„  +  a.y.)) 
i  i  11  Z  Z  J  j  4  4 


+  k2(v2  -  (b1y1  +  b2y2  +  b^  +  b^))  +  (o )  +  (q ) 


+  A?(0)  +  Ag(0)  +  Ag(o)  +  A10(0)  +  ^11(0)  +  *12(0) 


+  A13  (-  tT  (y7-y5))  +  A14  (-  T 
b  e 


+  A15  T  (G2G3  Vry7))  +  (viai  +  °ibi> 
g  1=1 


The  necessary  and  sufficient  conditons  for  minimum  by 
Pontryagin's  Minimum  Principle  are  [ t ] : 


■ 
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Y  =  - 


9H 

9A, 


:  9h 

V  w 


(5.2.5) 


:  9h 

V  3A 


:  _  9h 

A3  3B 


_9H 

9V 


=  0 


i_H  >  0 

2 
9  V 


where 


Y  =  [yl  y2  y3  y4  y5  ^6  y7 1 


A1  [X1  X2  X3  X4  A13  A14  A15] 


A2  [A5  A6  A7  V 


A3  [A9  A10  A11  A12^ 


and  the  boundary  term  is  minimized  with  respect  to 


y.(t£),  a  <t  ),  a  (t  )  b1(tQ) ,  b.(tf) 


which  yields 


■ 
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^i(tf)  =  0 


(1=1,2,3,4,13,14,15) 


(5.2.6) 


A  (t  )  =  A  (t  )  =0  (j=5, 6, 7, 8, 9, 10, 11, 12) 
Jo  J  I 


The  optimizing  equations  (5.2.5)  are  arranged  and  given 
specifically  as: 


h  =  X4y2  +  y2y4 


y2  =  S  !y5  ‘  S4V’l 


S5X3yl  "  Kdy2  “  S5Xly3  "  S4Xly4  S4yly4 


-  S  y.y0  +  K [u (t-t  )  -  u(t-x)]] 
CL  1  3  o 


y3  “  y6  '  A  y3  +  C  yl 


y4  ■  -  -  xiy2 


y5  =  J~  [y7  '  y5] 
b 


y6  =  7“  [crv2  '  y6] 

e 


y7  =  j-  [g2g3Vi  '  y7] 

g 

A  A 

\  -  2“yx  +  \y2  +  F  S5y3  +  F  S4y4  '  Vl  "  k2bl 


A  A 

+  F  s4x4  +  F  S5*3 


' 
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*  /  2 
A2  =  2a2y2  +  A4yl  "  Aly4  "  kla2  "  k2b2  +  X1A4  "  X4A1  +  M~  Kd 


A2  A2  s 

A3  =  2a3y3  +  M~  S5yl  ”  kla3  "  k2b3  +  A  A3  +  M~"  S5X1 


A  A 

A,  =  2ay.  -  A  y  +  —  S,y.  -  k  a  -  k  b  -  C  A  +  — -  S.x^ 
4  4  12  M  4  1  14  24  3  M41 


A5  =  2vlal  "  klyl 


A6  =  2v2a2  "  kly2 


X7  =  2V3a3  -  kly3 


A0  =  2v,  a,  -  k-y. 
8  4  4  1  4 


(5.2.7) 


A9  =  2a1b1  -  k2y1 


A10  =  2°2b2  -  k2y2 


A11  =  2°3b3  '  k2y3 


X12  =  2°4b4  -  k2y4 


A13  A2 
A13  “  2yly5  +  T,.  ”  M 


A 


14 


A14  =  2y2y6  +  '  X3 


- 
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\ 


A15  =  2Y3y7  +  T 


15 


13 


g 


.  G2G3  X15 
V1  2  6,  T 


1  g 


2  6, 


c  xi/ 

r  14 


v2  2  B0T 


2  e 


2  6, 


Equations  (5.2.4),  (5.2.6)  and  (5.2.7)  constitute  a  boundary 
value  problem  for  the  system  and  these  are  to  be  solved  for  the 


feedback  parameters  of  the  system. 


5s 


- 


CHAPTER  VI 


COMPUTATIONS  OF  THE  SYSTEM  EQUATIONS 

6 . 1  Computations  for  Turbo-Generator 

Equations  (5.2.7)  and  the  feedback  parameters  are  solved 

by  the  following  method. 

Equation  (5.2.3)  can  be  reduced  to 


Y  =  F(y,v,t)  ,  Y(tQ)  =  0 


(6.1) 


V  =  Z(a,b,y)  (6.2) 


where 


y  =  (y1,y2,y3,y4,y5»y6»y7) 


v  =  (vrv2) 


a  —  (a.^ , a2 » a^  , a^) 


b  =  (b1,b2,b3,b4) 


The  feedback  parameters  (a,b)  of  the  system  are  two  constant 
vectors  to  be  determined  such  that  the  cost  function  equation 

(5.2.2) 


J  -  J  £  [a(y£  +  y 4)  +  “2y2  +  a3y3  +  Yly5  +  Y2y6  +  Y3y7 
t 


+  B  v*  +  +  l  (V^  +  a.b^)]dt 

i=l 


=  /  f  L(y,  v,  a,  b,  t) 


dt 


(6.3) 
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is  minimized  subject  to  Y  =  F(y,v,t). 

This  feedback  control  problem  can  be  solved  by  a  simple 
reformulation  of  the  above  problem.  The  reformulation  is  carried 
out  by  substitution  of  equation  (6.2)  into  equation  (6.1)  and 
equation  (6.3),  which  gives 

tf 

J(a,b)  =  tJ  L[y,Z(a,b,y) ,  a,b,t]dt 
o 

The  fact  that  Y  are  the  known  functions  results  in 

tf 

J(a,b)  =  /  L(y ,a,b,t)dt 

t 

o 

Y  =  F(y,a,b,t) 

Now,  this  feedback  control  problem  becomes  a  problem  of 
obtaining  the  values  of  parameters  (a,b),  which  minimize  the 
cost  functional  subject  to  the  associated  differential  constraints. 
Thus  we  have  formulated  a  problem  of  parameter  optimization 


* 
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subject  to  Y  =  F(y,a,b,t).  The  feedback  parameters  should  satisfy 
the  two  point  nonlinear  boundary  value  problem  according  to 
equation  (5.2.7)  . 

To  solve  this  nonlinear  two  point  boundary  value  problem,  one 
has  to  resort  to  iterative  numerical  methods  (Section  2.4)  for  obtain¬ 
ing  the  optimal  feedback  parameters  and  trajectories  on  the  digital 
computer . 

The  first  variation  (Section  2.2)  of  the  functional  is 
expressed  as: 


o 


=  /  f  [hV  6a1  +  (H1)1  6b1] dt 

^  a  b 

o 

The  gradient  algorithm  for  the  cost  functional  with  respect  to 
the  feedback  parameters  are  given  by  [3] 


i+1 

a 


a  g 


,  i+1  ,i  ,  i 
b  =  b  -ah 


where 


g1  -  ^  -  H1 

3a 


l 

i  a 


h1  -  ^ 

3b1  b 


Hence,  if  the  change  in  a  ,  b  is  selected  as 
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x  i  i+1  i 
oa  =  a  -  a 


a 


5H~ 

3a‘ 


6b 


i 


3H 
a  — 

3b 


i 


i 


with  a  >  0  ,  then 


o 


If  these  nominal  feedback  parameters  satisfy  —  (y1,a1,b1, A1, t)  -0 » 

d  3. 

then  a1,b1  are  the  desired  optimal  feedback  parameters  a*,b*. 

For  a  >  0,  since  the  increment  in  is  governed  essentially 
by  the  first  variation  6j^,  this  implies 


_  r  i  -  i,  _  r  i+1  ,  i+1,  ,r  /N 

J-^ [a  ,b  ]  >  [a  ,b  ]  (6.4) 


Hence  from  the  equation  (6.4),  if  the  first  variation  ||6J^||<e  is 
satisfied  the  feedback  parameters  will  eventually  converge  to  the 
optimal  feedback  parameters. 

The  gradient  algorithm  for  the  optimum  are  then  given  by 


a1+1  =  a1  -  ag1  (j=l,2,3,4) 

J  J  J 

bi+1  =  bf  -  ahi  (j=l,2,3,4) 

J  J  J 


(6.5) 
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Also,  the  gradient  equations  are  obtained  as 
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Sj 


=  2v 


-  Vj 


(j-1,2,3,4) 


hj  “  2°jbj  "  k2yj  0  =  1, 2, 3 ,4) 

where 

g2g3 

kl  ■  T  A15  ■  2ei(alyl  +  a2y2  +  a3y3  +  a4y4> 

c 

k2  =  Y~  X14  “  2e2('blyl  +  b2y2  +  b3y3  +  b4y4^ 
e 

The  steps  to  determine  the  optimal  feedback 

parameters  a*,b*  on  the  digital  computer  using  the  steepest  descent 
method  are: 

Step  1: 

Consider  the  parameters  a1,b1  as  being  a  piecewise-constant 

during  t  <  t  <  t_.  Assume  the  initial  values  a1  =  a  ,  b1  =  b  . 

o  —  —  f  o  o 

Step  2: 

Using  a^’jb1,  integrate  the  state  equations  (5.2.3)  forwards  from 

t  to  t„  with  the  initial  conditions  y1(t  )  =  y  and  store  the 
of  o  ■'o 

resulting  state  trajectory  yX(t) . 

Step  3: 

Using  the  stored  value  of  yX(t) ,  integrate  the  costate  equations 

(5.2.7)  backwards  from  t.  to  t  with  the  final  condition 

f  o 

•  i  i  i  i 

Xx(t  )  =  X  and  compute  g  =  H  ,  h  =  H,  and  store  these  functions, 
f  f  a  d 

Step  4: 

Compute  the  first  variation 


J1[ai,b1]  -  JXU 


i+1 


b1+h 


_ 
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Step  5: 

If  the  stopping  criterion  | | SJ^ | |<e  is  satisfied,  terminate  the 
iterative  procedure.  Then  these  a"\b^  are  the  required  optimal 

•  |  *i  •  |  *i 

feedback  parameters.  Otherwise  generate  the  new  a1  jb1  according 
to  equations  (6.5),  and  update  a^,b^  by  a^"^,b^"^  and  return  to 
Step  2. 

6 . 2  Discussion  of  Results 

This  thesis  is  aimed  at  obtaining  a  desired  performance  of  the 

turbo-generator  during  the  period  following  a  distrubance  by  means  of 

an  optimum  variation  in  the  control  inputs.  The  objectives  of 

optimization  are  to  improve  the  stability  and  to  damp  out  any  oscillations. 

Both  these  objectives  can  be  achieved  by  using  a  suitable  function 

s  2 

of  the  rotor  angle  of  the  form  a (6  -  6  )  ,  where  a  is  a  constant 
penalty  factor.  As  shown  in  Chapter  5,  the  usual  quadratic  cost 
function  in  the  machine  variable  6  can  effectively  be  replaced  by  a 
new  cost  function  which  is  quadratic  in  the  variables  sin6  and  cos5. 

In  order  that  the  rotor  angle  converges  to  a  final  steady  value 
without  large  steady  state  errors,  the  cost  function  is  taken  as  a 
quadratic  function  of  the  deviation  of  the  rotor  angle  from  the  final 
steady  state  value.  A  function  of  the  field  flux  linkage  deviation 
a3(x3  -  x ^ )  ,  where  is  a  constant  penalty  factor,  is  also  included 
to  counteract  the  tendency  of  the  rotor  flux  to  decay  during  the 
disturbance  period.  This  helps  to  reduce  steady  state  errors  in 
the  terminal  voltage  of  the  machine  at  the  end  of  the  optimization 
period . 

6.2.1  Choice  of  Penalty  Factors 


Penalty  factors  are  used  in  the  cost  function  in  order  to 


' 

s 
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specify  the  relative  importance  of  the  various  terms  which  are 
minimized.  In  this  thesis,  penalty  factors  are  chosen  so  that 
the  major  proportion  of  the  reduction  in  the  value  of  the  cost 
function  in  each  iteration  is  obtained  in  the  term  containing  6. 
a  and  are  chosen  so  as  to  alter  the  cost  function,  and  to 
compensate  for  the  change  due  to  the  transformation  as  used  in 
Chapter  4.  Other  penalty  factors  are  chosen  from  [2]. 

6.2.2  Performance  of  Turbo-Generator 

The  machine  is  disturbed  from  the  steady  state  by  a  torque  pulse 
applied  directly  to  the  shaft  and  the  generator  is  controlled 
through  a  linear  feedback  of  the  state  variables.  It  is  assumed  that 
the  disturbance  is  modelled  by  introducing  a  torque  pulse  in  the  right 
hand  side  of  the  torque  equation  of  magnitude  K  which  is  increased 

g 

in  steps  of  5%  of  steady  state  value  and  a  duration  t  (0.5  sec). 

Hence  in  these  results,  the  values  of  the  variables  at  t  =  t  are 

o 

the  steady  state  values. 

The  optimum  controls  and  the  optimum  swing  curves  are  shown  in 
Fig.  6.1,  Fig.  6.2,  Fig.  6.3,  Fig.  6.5,  Fig.  6.6,  Fig.  6.8,  Fig.  6.9 
and  Fig.  6.14.  As  shown  in  Fig.  6.1,  Fig.  6.2  and  Fig.  6.3,  the 
variation  of  the  rotor  angle  is  increased  as  the  disturbance  strength 
K  is  increased. 

When  the  disturbance  is  at  K  =  0.4,  the  control  signal  to  the 
speeder  gear  drops  near  zero  value  requiring  the  valves  to  close 

in  Fig.  6.9.  This  is  followed  by  a  steady  signal 
at  the  frequency  of  the  rotor  oscillations  5  the  nature  of  the 
variations  depending  on  whether  the  rotor  is  accelerating  or  de¬ 
celerating.  The  exciter  voltage  control  at  K  =  0.4  is  shown  in  Fig. 


« 
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6.8.  It  is  shown  that  a  suitable  control  at  K  =  0.4  can  be  found 

to  improve  the  dynamic  performance  of  a  large  turbo-generator,  even 

in  the  presence  of  practical  constraints  on  the  control  system. 

The  swing  curve  A  in  Fig.  6.3  shows  instability  at  K  =  2.1.  It 

is  of  concern  where  a  large  disturbance  occurs  without  any  controls. 

However,  the  swing  curve  B  in  Fig.  6.3  shows  stability  with  optimum 

feedback  controls.  This  system  is  stable  since  it  returns  to  a 

state  of  equilibrium  when  disturbed  from  a  state  of  equilibrium. 

Fig.  6.4  shows  a  phase  plane  plot  corresponding  to  Fig.  6.3. 

The  plot  A  is  unstable  and  the  torque  angle  6  passes  unstable 

equilibrium  point  6^  without  the  controls.  However  the  motion  B 

is  governed  by  optimum  feedback  controls  and  the  accelerating  toraue 

(K  =  2.1).  The  machine  swings  to  accelerate,  reaching  its  maximum 

position  6  and  starts  to  decelerate,  then  returning  towards  6  . 

m  f 

The  oscillations  will  finally  stop  the  rotor  at  6  with  damping  after 

s 

the  disturbance  period  (0.5  sec.).  Fig.  6.7  also  shows  a  phase 
plane  plot  at  K  =  0.4. 

In  conclusion,  the  optimum  feedback  controls  are  used  for  the 

improvement  of  turbo-generator  stability  after  a  disturbance. 

6.2.3  Effectiveness  of  Optimum  Feedback  Controls 

The  parameters  of  the  system  are  the  same  as  [2]  and  the 

disturbance  period  is  0.5  sec.  The  exciter  parameters  used  are 

Te  =  0.2  sec.  and  C  =1.0.  An  exciter  system  which  consists  of  an 

r 

exciter  or  transformer  with  thyristors  is  assumed,  where  a  fast 
rate  of  response  can  be  obtained.  The  gain  can  be  taken  as  unity 
without  loss  of  generality  since  open  loop  control  is  assumed  and 
per  unit  values  are  used  throughout. 
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Fig.  6.1  shows  the  swing  curves  at  K  =  0.4  and  Fig.  6.8, 

Fig.  6.9  show  the  optimum  feedback  controls  required.  Fig.  6.10  and 
Fig.  6.11  show  the  corresponding  variations  of  the  field  voltage 
and  the  mechanical  torque.  Fig.  6.13  shows  the  corresponding 
variation  of  the  excitation  voltage. 

The  high  speed  voltage  regulator  and  exciter  have  two  major 
effects.  The  first  is  the  increase  in  restoring  synchronizing  forces 
made  possible  through  the  forcing  of  excitation  and  internal  machine 
fluxes.  This  improves  transient  stability.  The  second  is  the 
deterioration  of  machine  damping,  resulting  in  dynamic  instability. 
However,  it  is  possible  to  produce  some  effective  control,  through 
the  variation  of  excitation  as  shown  in  Fig.  6.13,  Fig.  6.10  and 
Fig.  6.8. 

Fig.  6.9  shows  the  governor  action  at  K  =  0.4  and  Fig.  6.11 
shows  the  corresponding  variation  of  mechanical  torque.  In  this 
investigation,  the  governor  and  relay  combined  time  constant  is  taken 
as  0.2  sec.  and  the  turbine  time  constant  as  0.49  sec.  In  actual 
practice,  the  governor  and  relay  time  constants  could  be  smaller 
depending  on  the  type  of  hydraulic  equipment;  and,  in  fact,  the  turbine 
time  constant  is  only  representative  of  the  high  pressure  cylinder 
delay  in  a  compound  turbine.  This  is  effective  in  preventing 
excessive  speed  variations  during  the  disturbance,  and  some  effects 
are  possible  with  fast  response  as  shown  in  Fig.  6.9,  Fig.  6.11  and 
Fig.  6.12. 

Fig.  6.5  and  Fig.  6.6  show  the  trajectory  of  rotor  angle  after 
0.5  sec.  The  swing  curves  B  in  Fig.  6.5  and  Fig.  6.6  are  with 
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optimum  feedback  controls.  Fig.  6.15  and  Fig.  6.16  show  the  trajectory 
of  optimum  feedback  parameters.  The  optimum  feedback  parameters 
are  time-independent.  However,  it  is  shown  that  they  depend  on  the 
disturbance  strength  K. 

The  curve  A  in  Fig.  6.14  shows  the  swing  curve  at  K  =  0.5 
with  optimum  feedback  parameters  a_^,  b.  found  at  K  =  0.5  and  the 
curve  B  shows  the  swing  curve  at  K  =  0.5  with  optimum  feedback 
parameters  a.,  b.  found  at  K  =  0.2.  The  cost  function  on  the  curve 
A  is  J  (opt)  =  0.50186658,  but  the  cost  function  on  the  curve  B  is 
J  (cost)  =  0.51556063.  Replacing  the  parameters  (K  =  0.5)  by  the 
parameters  (K  =  0.2)  at  K  =  0.5,  the  value  of  the  cost  function 
becomes  higher. 


ROTOR-RNGLE(DEG) 
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FIG.  6.1  ROTOR-RNGLE  VS  TIME  (  RT  K=  0.4  ) 


ROTOR-RNGLE(DEG) 
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FIG.  6.2  ROTOR-RNGLE  VS  TIME  (  RT  K=  1.0  ) 


ROT  OR- ANGLE(DEG) 
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FIG.  6.3  ROTOR-ANGLE  VS  TINE  (  AT  K=  2.1  ) 


RNGULRR-VELOCITY  (rad/sec) 
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FIG.  6.4  PHRSE  PLANE  PLOT  AT 


K=  2.1 


ROTQR-FINGLE(DEG) 
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TIME  (SECOND) 

FIG,  6.5  ROTOR-ANGLE  VS  TIME  (  AT  K=  0,1  ) 


ROT  OR-RNGLE(DEG) 
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FIG.  e.e  ROTOR-ANGLE  VS  TINE  (  AT  K=  0.4  ) 


ANGULAR-VELOCITY  (rad/sec) 
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ROTOR-flNGLE(DEG) 


FIG.  6-7  PHASE  PLANE  PLOT  AT 


K=  0.4 
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D 


FIG.  6.8  OPTIMUM  CONTROL  TO  TURBO~GENERRTOR  MODEL 


FIG.  6.9  OPTIMUM  CONTROL  TO  TURBO-GENERATOR  MODEL 

(K=0.4) 
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o 


TIME  (SECOND) 

FIG.  6.10  VARIATION  OF  OPTIMUM  FIELD  VOLTAGE  (k«o.4) 


o 


FIG.  6.11  VARIATION  OF  OPTIMUM  MECHANICAL  TORQUE 

(K=0 .4) 
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l/D 


TIME  (SECOND) 


FIG.  6.12  VARIATION  OF  OPTIflUfl  ANGULAR  VELOCITY 

(K=0 .4) 
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o 


FIG.  6-13  VRRIRTION  OF  OPTIMUFI  EXCITRTION  VOLTRGE 

(  K-  0 .  A  ) 


ROT  OR-flNGLE(DEG) 
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TIME  (SECOND) 

FIG.  6.i4  ROTOR-ANGLE  VS  TINE  (  RT  K=  0.5  ) 


optimal  PARAMETERS  (AI) 
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FIG.  6.15  OPTinflL  PRRRMETERS  VS  DISTURBRNCE 


optimal  PARAMETERS  (BI) 

.120  -0-060  0.000  0.060  0.120 
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FIG.  6.16  OPTIMAL  PARAMETERS  VS  DISTURBANCE 


CHAPTER  VII 


CONCLUSIONS 

The  optimizing  equations  for  a  turbo-generator,  which  is 
connected  to  an  infinite  bus  and  controlled  through  a  linear 
feedback  of  the  state  variables,  are  obtained.  It  should  be 
noted  that  the  transfer  functions  of  the  governor  and  the  exciter 
are  included. 

The  optimal  feedback  parameters  are  obtained  by  applying  the 

gradient  descent  method  to  the  two  point  nonlinear  boundary 

value  problem.  The  values  to  be  obtained  for  these  parameters 

depend  on  the  strength  K  and  the  duration  t  of  the  disturbance 

since  the  model  is  nonlinear  contrary  to  the  usual  feedback  control 

of  a  linear  model.  It  is  assumed  that  the  strength  K  of  torque 

pulse  is  increased  by  5%  of  steady  state  value  and  the  duration  t 

of  torque  pulse  is  0.5  sec.  in  the  program. 

Using  the  optimal  feedback  parameters  the  nonlinear  state 

equations  on  the  feedback  are  integrated  by  using  Runge-Kutta 

method  fowards  from  t  =0  to  t  =0.5  sec.  and  then  the  nonlinear  costate 

o  f 

equations  are  integrated  by  using  the  same  method  backwards  from 

t  =0.5  to  t  =0  sec.  The  value  of  the  cost  functional  is  also 
f  o 

obtained  by  Simpson’s  method. 

In  this  research,  the  optimum  is  reached  in  less  than  11 
iterations  and  within  a  computing  time  of  1  minute  on  the  AMDAHL  470/V7. 

The  extension  of  the  method  to  a  multimachine  system  is 
feasible  by  optimizing  each  machine  at  every  time  step.  Although 
this  increases  computation  time,  the  complexity  of  the  method  is  not 
changed.  The  sensitivity  of  the  system  to  parameter  variations  and 
the  effects  of  damper  windings  are  also  of  concern. 
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A.  Steady  State  Values  of  Machine  Parameters  Machine 

Generator  rating:  37.5  MVA  or  30MW  at  0.8  PF. 

Bus  voltage  (line-line)  :  1  p.u. 

No.  of  poles  :  2 

R.P.M.  =  3000 


f  =  50  Hz 


x  =  0.1265  p.u. 
t  r 

e 

on  100  MVA 


x  =  0.354  p.u.  on  100  MVA 
t 

r 

ratio  =  11.8/141.6 


Base : 

One  changes  from  100  MVA  base  to  37.5. 


x  =  0.354  x 

t 

37.5 

100 

p.u. 

=  0.1328 

p  .u . 

r 

x’  =  0.1265 

t 

37.5 
X  100 

p.u. 

=  0.0475 

p  .u . 

e 


The  diagram  becomes: 


Constants : 


infinite  bus 
141.6  KV 


V  =  1/0°  p.u. 


;* 
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Vb  is  the  line  to  ground  phasor  voltage 

11.8 

vb  =  ft  =  6.83  KV 


2_H  _  2x6.63  _  2 

co  2tt  x  50  ~  0 . 04225  sec  /rad 


Kd  =  0.02535 


Rf  =  0.00107  p.u 


x  =  0.14  p.u. 
a 


md 


=  1.86  p.u. 


x  =  1.86  p.u. 
mq 


x^d  =  2.00  p.u. 


x  =  x ' 


't  +  x^  =  0.133  +  0.048  =  0.181  p.u 


r  e 


V,  -  E  =  1.0  p.u. 
b  rms  r 


S,  = 


-vl  (2. 18 1-2. 18 lx  (1--y6---) 
b _ 2 _ 

2 

2.181(2.181-  (1,v6)--) 


=  -  1.76 


2 


■- 
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Sr.  = 


Vmd 


Jl 


( 


Lfd  x,+x 
d  e 


md 


‘fd 


0.707  x  ( 


r)  =  1.46 


2.181  - 


1.86 


A  = 


R,-  co  (x,  +  x  ) 
f  o  d  e 

x2h 

x  (x  +  x - — ) 

fd  d  e  x  ' 
id 


=  Q « 00107  X  50  X  2tt(2 ,181) 

2 

2.0  x(2  .181  -  1^-6  ) 


=  0.812 


C  = 


/2  V,  R_  to  x  , 
b  f  o  md 

xfd(x  +  X  - 

fd  d  e  x,, 

f  d 


v^2"  x  0.00107  x  2 7T  x  50  x  1.86 

2 

2.0  x  (2.181  -  1^6  ) 


=  0.98 


Since  PF  =  0.8,  I  =  0.8  -j0.6 

a  J 


Eq  =  1.00  +  j  0.0  +  j  (2 . 181)  (0.8  -  j  0. 

=  2.31  +  j  1.75 

=  2.89  /37°.  1 

Therefore  6  =  37°. 1 


At  the  steady  state  value  of  flux, 


d6 

dt 


0 


d(w  \p  ) 
o  s 

dt 


0 


co  V r  C  cos6 

i  /  \  oY  f 

“oVs>  = - x — — 

=  100tt(2 .35x10~3)  +  0.98cos37°  .1 

0.812 


=  1.874. 


From  the  torque  equation, 
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A 

M  (s)  =  -r—  sin2  6  +  go  ip  (s)  sin 6 

t  Z  S  -)  O  I  s 


-1.76 

2 


sin(2x37 ° .1)  +  1.46  x  1.874  x  sin37°.l 


=  0.8 


Vf(s) 


E  R. 
o  f 


x 


md 


yfj  x  2.89  x  0.00107 

1.86 


=  2.35  x  10 


p  .u . 


Let  go  ib- 
of 


GO  Vr 

o  f 


* 


then  the  equations  are: 


— - —  sin6  cos6 

M  M 


sin  6 


aa 

M  dt 


+  C  cosS 


where 


M  =  0.04225  p.u. 

K, 

=0.02535/0.04225  =  0.6  p.u. 

M 

s  / 

— —  =  -1.76/0.04225  =  -41.65  p.u. 
M 


•• 


4 
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=  1.46/0.04225  =  34.5  p.u. 

A  =  0.812  p .u . 

C  =  0.98  p.u. 

6  (s)  =  37°. 1 


x^(s)  =  0.0 
u|(s)  =  0.8 

u^(s)  =  2.35xl0"3  x  IOOtt  =  0.738 

u, (s)  =  Y  (s)  =  0.424 

1  o 

u_  (s )  =  co  V  (s)  =  0.738 

2  o  ex 

Also  these  parameters  are  used  in  all  the  optimal  feedback  control 
studies. 

G2  =  1.33 

G3  =  1.42 

T  =0.2  sec. 
g 

T,  =  0.49  sec. 
b 

C  =  1.0 
r 

T  =0.2  sec. 
e 

For  the  cost  function  penalty  one  takes: 
a  =2.50 
a2  =1.00 


0.10 


99 


$  =  1.00 

e2  =  1.00 

y  =  o.oi 

Y2  =  0.01 

Y3  =  0.01 

v±  =  0.001 

v2  =  0.001 

v3  =  0.001 

v  =  0.001 

4 

o  =  0.001 

a 3  =  0.001 

a3  =  0.001 


0.001 
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C.  Expression  of  Terminal  Voltage  in  Terms  of  State  Variables 


Fig.  C.l  : 

u 

1.0  +  j  0.0  (p 

rms 

V 

/2  U 

e 

rms 

V  ,  = 

/l  U  sin  6 

ed 

rms 

V 

/2  U  cos  6 

eq 

rms 

Under 

the  assumptions 

:  Hi  0, 
a 

we  have 

V  =  V 

,  -  x  i 

d 

ed  e  q 

V  =  V 

-  x  i  , 

q 

eq  e  d 

a, 


0.  iq  £  0,  re  Hi  0 


(C.l) 


V 


d 


a)  ip 

o  q 


(C .  2) 


V  =  -  u)  iK 
q  o  d 
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to  ^  ,  =  x  i  +  x  i 
o  d  d  d  ad  f 


co  =  x  i 
o  q  q  q 


(C .  3) 


Substituting  (C.3)  into  (C.2),  it  becomes 


V  -  x  i 
d  q  q 


V  =  “  x  ,  i  -  x  i 
q  d  d  ad  f 


(C .  4) 


where  is  the  terminal  d-axis  voltage  and  is  the  terminal  q- 

axis  voltage.  Using  co  ip-  =  x  .i.  +  xriri  we  have 

of  ad  d  it 


V  =  -  x  ,i  ,  = 

q 


ad  /  ,  .  \ 

J-J  - —  (co  ^  -  X  1  ) 

d  d  x__  of  ad  d 
f  d 


and  hence  (C.4)  becomes 


V  =  x  i 
d  q  q 


x  ,  x 

V  =  - x  )  i, - ~  <*>  ip 

<>  xfd  d  d  xfd  0  f 


Solving  for  i^  and  i^,  we  obtain 


i  =  (v  +  to  ip  )  /  (—  -  x  ) 
d  q  x£,  o  f  x_  d 


fd 


fd 


V 

.  c 

l  =  — 

q  x. 


(C  .5) 
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Substituting  (C.5)  into  (C.l),  we  obtain: 


/2  U  x 
TT  rms  q 

V,  =  - — - ^  sin  6 

d  x  +  x 
q  e 


x  -  x.  xr  , 

TT  _  A7  TT  ad  d  fd 

V  =  v2  U  — n-  -  - - - - —  cos  6 

q  rms  2  ,  . 

x  -  (x  -+x  )  x 

ad  d  e  fd 


+ 


x  x  _  j 

e  f  d 


x  -  (x  -he  )  x 
ad  d  e  fd 


Vf 


Data:  x  =  2.0,  x  =  2.0,  x  =  0.181,  x  ,  =  1.86,  x£J  =  2.0 
d  q  e  ad  fd 

6  =  37.1°,  a )  =  -  1.874 

o  f 

Solving  for  and  by  using  the  datas,  we  obtain: 


Vd  =  0.782  (p.u.) 
V  =  1.374  (p.u.) 

q 


Therefore 


J  +  v2 

V  =  — 2 - =  1.118  (p.u.) 

t  — 


(C .  6) 


/2 


M 
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Id 


Fig.  C.2  Generator  Phasor  Diagram 


Therefore 


V  =  E  -  j  I  x 
t  f  J  a  q 

=  E  +  j  i  x 

e  a  e 

=  (1.0  +  j  0.0)  +  j  (0.8  -  j  0.6)  X  0.181 


=  1.1086  +  j  0.1448 

=  1.118  /7 ♦ 44°  (C  . 7 ) 

Comparing  (C.6)  with  (C.7),  the  generator  terminal  voltages  are 
the  same.  Thus,  the  generator  terminal  voltage  depend  on 
co  ip  and  6. 
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D .  Results 

It  is  assumed  that  the  strength  K  of  torque  pulse  is  increased  by 
5%  of  steady  state  value  and  the  duration  t  of  torque  pulse  is 
0.5  sec.  in  the  program. 

A.,  B.,  Y.,  P.  and  V.  in  the  results  are  corresponding  to  a.,  b., 
1111  l  ii 

y..  A.,  and  v_^  in  the  body  of  the  thesis  respectively. 


"  .  a 
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REALISTIC  FEEDBACK  CONTROL  OF  TURBOGENERATORS; 

THE  ALTERNATER  IS  CONTROLLED  THROUGH  A  LINEAR  FEEDBACK  OF 
THE  STATE  VARIABLES. 

THE  FEEDBACK  PARAMETERS  ARE  OBTAINED  BY  SOLVING  A  TWO  POINT 
NONLINEAR  BOUNDARY  VALUE  PROBLEM. 

THE  VALUES  TO  BE  OBTAINED  FOR  THESE  PARAMETERS  DEPEND  ON 

THE  STRENGTH  AND  DURATION  OF  THE  DISTURBANCE  SINCE  THE  MODEL  IS 

NONLINEAR  CONTRARY  TO  THE  USUAL  FEEDBACK  CONTROL  OF  A  LINEAR  MODEL 

DIMENSION  TBL (41 ,25) ,Y(88) , DY ( 88 ) , A ( 4 ) ,B(4) , AA (  4 )  , BB(4) ,V(2) , 

JCOST ( 3 1 ) , PLOT ( 2 1 ,8) 

REAL  K,KK, JCOST, JOPT.K1 ,K2 
INTEGER  SW 

DATA  A,B,AA,BB/16*0.0/,JCOST/31*0.0/ 

DATA  A/0.01633753,-0.41949201 ,-0.00091682,-0.01575423/, 

*  B/- 0.004538 19,0. 11382747,0.0,0.00430573/ 

K=-0 . 05 
DO  99  11=1,11 
K=K+0 . 05 


SW=  1 

CALL  PARAMT ( A , B  ,  K  ) 

' PARAMT'  IS  THE  SUBROUTINE  WHICH  OBTAINS  THE  OPTIMAL  FEEDBACK 
PARAMETERS  BY  SOLVING  THE  TWO  POINT  NONLINEAR  BOUNDARY  VALUE  PROBLEM. 

DO  10  d=1 ,88 
Y  ( <J )  =0 . 0 
DY ( J ) =0 . 0 

INITIAL  Y , DY  WITH  ZERO. 

DO  15  J=1 ,25 
DO  15  1=1,41 
TBL ( I , d ) =0 . 0 
INITIAL  TBL  WITH  ZERO. 


N  =  22 
T  =  0 . 0 
T  F  =  0 . 5 
L  =  0 
M=  1 
dd  =  0 

N  IS  NO.  OF  EQ. 

T  IS  INITIAL  TIME. 

TF  IS  FINAL  TIME. 

LISA  PARAMETER  IN  "RUNTA" . 

M  IS  NO.  OF  INTERVAL, 
dd  IS  COUNTER. 

CALL  OPT IM ( Y , DY , A , B , SW , K ) 
SUBROUTINE  OF  OPTIMAL  EQUATIONS 
CALL  RUNTA(7,L,I,Y,DY,T,0.005) 
SUBROUTINE  OF  RUNGE-KUTTA  METHOD 
GO  TO  (20,21 ), I 
U  J  =  d  J+ 1 

IF (MOD ( Jd , 10) )  20,22,20 
M=M+1 

CALL  VVKK ( Y , A , B , V , K 1 ,K2) 

SOLVE  CONTROLS  VI ,V2. 

DO  31  d=1 ,7 


|, *J tmahar 
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'  TBL  <  M , J ) = Y ( d) 

DO  32  d  =  23 ,24 
)  TBL ( M , J ) = V ( J-22) 

JCOST  (M)=COST(Y,A,B,V) 

IF ( M- 1 1 )  20,64,64 

50  THE  NON-LINEAR  STATE  EQUAT IONS ( Y 1  - -- Y7 )  ON  THE  FEEDBACK  ARE  INTEGRATED 
BY  USING  RUNGE-KUTTA  METHOD 

FORWARDS  FROM  T=0.0  SEC  TO  TF=0.5  SEC  BETWEEN  20  LINE  AND  50  LINE. 

1  CALL  SIMP ( JCOST , JOPT ) 

J ( OPT )  IS  OBTAINED  BY  USING  THE  SIMPSON'S  METHOD. 

CALL  EIGENV ( Y ( 8 ) ,Y(9),Y(11),EI,K1,K2) 

FIND  THE  SMALLEST  EIGENVALUE  OF  '  L'  . 

TBL ( 1 1 , 25 ) =EI 
d  d  =  0 

DO  65  J= 1 ,88 
Y( J ) =0 . 0 
5  DY ( d ) =0 . 0 

SW=2 
N=22 
T  =  0 . 5 
T  F  =  0 . 0 
L  =  0 
M=  1 1 

IF  SW= 1 , I T  IS  THE  FIRST  STEP  FOR  Y1---Y7. 

IF  SW=2 , I T  IS  THE  SECOND  STEP  FOR  PI P 1 5 . 

IF  SW=3 , 1 T  IS  THE  LAST  STEP  FROM  T  =  0.5  SEC  TO  TF  =  2.0  SEC. 

N  IS  NO.  OF  EQ. 

T  IS  THE  INITIAL  TIME. 

TF  IS  THE  FINAL  TIME. 

LISA  PARAMETER  IN  "RUNTA"  . 

M  IS  NO.  OF  INTERVAL. 


5  M=M- 1 

DO  68  J=1  ,7 
DY ( d ) =0 . 0 
5  Y(d)=TBL(M,d) 

)  CALL  OPTIM(Y,DY,A,B,SW,K) 

CALL  RUNTA (22,L,I ,Y,DY,T,-0.005) 
GO  TO  (70,71),I 
1  dd=dd+ 1 

I F ( MOD ( d d , 10) )  70,72,70 
l  dd  =  0 

DO  81  d=8 , 22 
I  TBL  ( M , J ) =  Y  (  d  ) 

IF ( M- 1 )  85,85,66 


THE  NON-LINEAR  COSTATE  EQUAT IONS ( P 1 - P15)  ARE  INTEGRATED  BY 

USING  THE  SAME  METHOD  BACKWARDS  FROM  T=0.5  SEC  TO  TF=0.0  SEC 
BETWEEN  66  LINE  AND  81  LINE. 

>  DO  86  d= 1 2 ,  1 9 

>  TBL ( 1 , d ) =0 . 0 


SW=3 


M 


KK  =  0 . 0 

IF ( SW . EQ . 3 )  GO  TO  94 
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DO  87  1=1,22 
Y ( I )  =0 . 0 
DY ( I ) =0 . 0 

DO  89  1=1,4 
Y ( I ) =TBL ( 11,1) 

AA( I ) =0 . 0 

BB( I ) =0 . 0 

T  =  0 . 5 
T  F  =  2 . 0 
M=0 
J=  1 1 
J  J  =  0 
L  =  0 

CALL  OPTIM(Y,DY,AA,BB,SW,KK) 
CALL  RUNTA(22,L,I,Y,DY,T,0.005) 
GO  TO  (90,91  )  ,  I 
J J  =  JJ+ 1 

IF ( MOD ( JJ , 10) )  90,911,90 


J  =  J+1 

DO  92  1=1,22 
TBL ( J , I ) =Y ( I ) 

CALL  VVKK ( Y , AA , BB , V , K 1 ,K2) 

DO  93  1=23,24 
TBL ( J , I ) = V ( I -22 ) 

IF ( 0-41 )  90,94,94 

T  IS  THE  INITIAL  TIME  TO  START  AFTER  THE  DISTURBANCE. 

TF  IS  THE  FINAL  TIME. 

M  IS  NO.  OF  INITIAL. 

J  IS  NO.  OF  TBL-ROW. 

FIND  OUT  WHAT  THE  DISTURBANCE  IS  AFTER  0.5  SEC 
BETWEEN  85  LINE  AND  93  LINE. 

CALL  WRITER(K,A,B, JOPT , TBL  ) 

CONTINUE 

STOP 

END 

NUMERICAL  SOLUTIONS  FOR  OPTIMAL  PARAMETERS 
USING  GRADIENT  TECHNIQUES; 

SUBROUTINE  PARAMT ( A , B , K ) 

DIMENSION  Y(88),DY(88),A(4),B(4),AK(4),BK(4),G(4),H(4),GG(4), 

*  HH  ( 4 )  ,TBL(  1 1 ,22)  ,  JCOST  (  3 1  ),  JOPT  (31  )  ,  J  J  ( 3 1  )  ,  V  ( 2  ) 
INTEGER  SW 

REAL  JCOST, K, JOPT, dd , K 1 , K2 , JK 1 , JK , JD 

DATA  X/1 . 0/ , TBL/242*0. 0/ , JCOST , JOPT/62*0. 0/ , EPSI/0. 0001/ 

*  , JJ/3 1*0.0/, JK/0.0/ 

Y  IS  THE  STATE  FUNCTIONS 
DY  IS  THE  DIFFENTIAL  FUNCTIONS 
A  IS  THE  FEEDBACK  PARAMETERS  A ( I )  A1,A2,A3,A4 
B  IS  THE  FEEDBACK  PARAMETERS  B(I)  B1,B2,B3,B4 


■ 
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AK  IS  A(I+1) 

BK  IS  B( 1  +  1  ) 

G  IS  DH/DA 

H  IS  DH/DB 

GG  IS  ( DH/DA  )  **2 

HH  IS  ( DH/DB ) **2 

JK  IS  THE  COST  FUNCTIONAL  J ( I ) 

JK 1  IS  THE  COST  FUNCTIONAL  J(I+1) 
JD  IS  d ( DELDA ) = J ( I ) - J ( 1+ 1 ) 

K  IS  THE  STRENGTH  OF  TORQUE  PULSE 
EPSI  IS  THE  EPSILON 
X  IS  THE  ALPHA 

DO  25  11=1 ,20 

DO  2  1=1,88 
Y ( I )  =0 . 0 
DY ( I ) =0 . 0 


T  =  0 . 0 
M=0 
d=1 
SW=1 
J  J  J  =  0 
L  =  0 

CALL  OPTIM(Y,DY,A,B,SW,K) 

CALL  RUNTA(7,L, I ,Y,DY,T,0.005) 

GO  TO  (5,6) , I 
JJJ=JJJ+1 

IF (MOD ( J J J , 10) )  5,66,5 
J  =  J+1 

DO  7  1=1,7 
TBL( J, I )=Y(  I ) 

CALL  VVKK ( Y , A , B , V , K 1 ,K2) 
dd( J)=COST( Y,A,B,V) 

IF ( J- 1 1 )  5,88,88 
CALL  SIMP ( J J , AREA ) 

UK  1 =AREA 

J  =  J-  1 

DO  9  1=1,7 
Y  ( I  )  =TBL ( J , I ) 

DY  ( I ) =0 . 0 

T  =  0 . 5 
M=0 
SW  =  2 
J  J  J  =  0 
L  =  0 

CALL  OPTIM(Y,DY,A,B,SW,K) 

CALL  RUNTA(22,L, I ,Y,DY,T, -0.005) 
GO  TO  ( 10,  1  1  )  ,  I 
jj j= j j j+ i 

IF  (MODI Odd  ,  10  )  )  10,111,10 

DO  13  1=12,19 
TEL ( J , I ) =DY  ( I  ) 

IF ( J-2 )  14,14,8 
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DO  16  1=12,15 
DO  15  0=1,11 
JCOST ( 0 ) =TBL (0,1 ) 

0OPT(0)=TBL(0,I)**2 
CALL  SIMP ( OCOST , AREA ) 

6(1-11 )=AREA/0.5 
CALL  SIMP ( OOPT .AREA ) 

GG ( 1-1 1 ) = AREA 

DO  18  1=16,19 
DO  17  0=1,11 
OCOST ( 0 ) =TBL ( 0  , 1 ) 

OOPT ( 0 ) =TBL ( 0 , I ) **2 
CALL  SIMP ( OCOST , AREA ) 

H ( I  - 15 ) =AREA/0 . 5 
CALL  SIMP ( OOPT , AREA ) 

HH ( 1-15) = AREA 

I F ( K . GT . 0 . 4  )  X  =  0.1 
IFIK.GT.  1  .3)  X  =  0 . 01 
IF ( K . GT . 2 . 0  )  X  =  0 . 01 

DO  20  1=1,4 
AK( I )=A( I ) -X*G(  I  ) 

BK( I )=B( I ) -X*H ( I  ) 

OD=OK-OK1 

DO  21  1=1,4 

IF  (  ABS ( G ( I ) ) .GT  .  E  PSI  .  AND.  ABS(  OD) .GT.EPSI  )  A ( I ) = AK (  I  ) 

IF(ABS(H( I ) ) .GT.EPSI .AND.ABS(OD)  .GT.EPSI  )  B ( I )  =  BK (  I  ) 

OK  =  OK  1 

IF ( ABS ( OD ) .LE.EPSI )  GO  TO  27 
CONTINUE 

I F ( 1 1 . GE  .  1 )  GO  TO  27 
DO  27  1=1,4 

IF ( SORT ( GG ( I ) ) .GT.EPSI )  GO  TO  1 
1F(SQRT(HH( I )) .GT.EPSI  )  GO  TO  1 
CONTINUE 

PRINT  30,11 
FORMAT (15) 

RETURN 
END 

OPTIMAL  EQUATIONS  Y1 - Y7,P1 - P15 

SUBROUTINE  OPT IM ( Y , DY , A , B , SW , K ) 

DIMENSION  Y ( 88 ) , DY ( 88 ) , A ( 4 ) ,B(4)  , V ( 2  ) 

INTEGER  SW 

REAL  M,KD,K,K1,K2,N1,N2,N3,N4 

DATA  M/0. 04225/, KD/O. 02535/, S4/-1 .76/.S5/1 .46/, AA/O . 8 1 2/ , C/0 . 98/ , 

*  X3S/1 .874/ ,U1S/0.8/,A1/2.5/,A2/ 1 .0/.A3/0. 1 / , A4/2 . 5/ , B 1 / 1 . 0/  , 

*  B2/1 . 0/ , X1S/0 .603/ .X4S/0.798/ ,C1/0. 01/ .C2/0.01 / ,C3/0. 01/  , 

*  G1/0. 00188/ ,G2/1 .33/.G3/1 .42/,TG/0.2/,TB/0.49/,TE/0.2/,CR/1  .0/ 
DATA  N1 , N2 , N3 , N4 , 0 1 , Q2 , Q3 , 04/8*0 .001/ 


CALL  VVKK ( Y , A , B , V , K 1 ,K2) 


' 
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VI ,Y2,Y3,Y4,Y5,Y6,Y7  OPTIMAL  EQUATIONS. 

IF ( SW . EQ . 2 )  GO  TO  10 

DY(  1  )=X4S*Y(2)+Y(2)*Y(4) 

DY(2)=1/M*(Y(5)-S4*X4S*Y( 1 )-S5*X3S*Y(  1 ) -KD*Y ( 2 ) -S5*X 1 S*Y ( 3 ) 

*  -S4*X 1 S*Y ( 4 ) -S4*Y ( 1 ) *Y ( 4 ) -S5*Y ( 1 ) *Y ( 3 ) +K*U  1 S ) 

DY ( 3 ) = Y  (  6 ) -AA*Y  <  3 ) +C* Y ( 4 ) 

DY ( 4 )  =  - Y ( 1 ) *Y ( 2 ) - Y ( 2 ) *X  1  S 
DY ( 5 )  =  1 / T B * ( Y ( 7 ) -Y ( 5 )  ) 

DY(6)  =  1/TE*(CR*V(2)-Y (6)  ) 

DY(7)=1/TG*(G2*G3*V( 1  ) - Y ( 7  )  ) 

P 1 , P2 , — . , P 1 4 , P 1 5  COSTATE  OPTIMAL  EQUATIONS 

IFISW.EQ. 1 )  GO  TO  20 

DY(8)=2.0*A1*Y(  1  )  +  Y( 1 1 )*Y(2)+Y(9)*S5*Y(3)/M+Y(9)*S4*Y(4)/M 

*  - K 1  *  A ( 1 ) -K2*B ( 1 )+Y(9)*S4*X4S/M+Y(9)*S5*X3S/M 
DY(9)=2.0*A2*Y(2)+Y( 1 1 )*Y( 1 ) - Y ( 8 ) *Y ( 4 ) -K 1 *A ( 2 ) -K2*B ( 2 ) +X 1 S*Y ( 1 1 ) 

*  -X4S*Y(8)+KD*Y(9)/M 

DY(10)=2.0*A3*Y(3)+Y(9)*S5*Y( 1 )/M-K1*A(3)-K2*B(3)+AA*Y( 10)+ 

*  S5*X 1 S*Y ( 9  ) /M 

DY ( 1 1)=2.0*A4*Y(4)-Y(8)*Y(2)+Y(9)*S4*Y( 1 )/M-K1*A(4)-K2*B(4) 

*  -C*Y( 10)+S4*X1S*Y(9)/M 
DY ( 12 ) =2 . 0*N 1 *A ( 1 ) - K 1  * Y (  1  ) 

DY ( 13)=2.0*N2*A(2)-K1*Y(2) 

DY( 14) =2. 0*N3*A( 3) -K1*Y ( 3) 

DY(15)=2.0*N4*A(4)-K1*Y(4) 

DY ( 16 ) =2 . 0*Q1 *B ( 1 ) -K2*Y(  1  ) 

DY(17)=2.0*Q2*B(2)-K2*Y(2) 

DY(18)=2.0*Q3*B(3)-K2*Y(3) 

DY ( 191=2. 0*Q4*B(4)-K2*Y(4) 

DY(20)=2.0*C1*Y(5)+Y(20)/TB-Y(9)/M 
DY(21)=2.0*C2*Y(6)+Y(21 )/TE-Y( 10) 
DY(22)=2.0*C3*Y(7)+Y(22)/TG-Y(20)/TB 
CONTINUE 

RETURN 

END 

CHECK  POSITIVE  DEFINITE  OF  '  L'  . 

SUBROUTINE  EIGENV ( P 1 , P2 , P4 , EIGEN , K 1 , K2 ) 

DIMENSION  A ( 17, 17) , WK ( 17) 

COMPLEX  W( 17) ,Z( 17, 17) ,ZN 

REAL  PI ,P2, P 4 , K 1 ,K2,S4,S5,M,N1 ,N2,N3,N4 

INTEGER  N, 1A, I  JOB , IZ, IER 

DATA  A1/2.5/.A2/1 .0/.A3/0. 1/,A4/2.5/,B1/1 .0/.B2/1 .0/.S4/-1  .76/, 

*  S5/1. 46/, M/0. 04225/, Cl, C2, 03/3*0.01/ 

DATA  N 1 , N2 , N3 , N4 , Q 1 ,Q2, 03,04/8*0 . 001/ 


I  A=  1 7 
IZ=  1 7 
N=  1 7 
I  JOB  =  0 

DO  10  d=1  ,  17 
DO  10  K=  1  ,  17 
A(  d , K )  =0 . 0 


A  ( 1  , 1 )  =  A 1 
A(  1  , 2 ) =  P 4 / 2 . 0 


. 

■ 
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A( 1 , 3 ) = ( P2*S5 ) / (2. 0*M ) 

A(  1 , 4 ) = ( P2*S4 ) / (2. 0*M ) 

A(  1 , 5) =-K1/2 . 0 
A( 1,9)  =  -K2/2 . 0 
A( 2 , 1 )  =  P 4 / 2 . 0 
A  ( 2 , 2 ) =A2 
A(2 , 4 ) =-P 1 /2 . 0 
A(2,6) =-K 1/2 . 0 
A ( 2 , 10)  =-K2/2 . 0 
A ( 3 , 1 )  =  (  P2*S5 ) / ( 2 . 0*M ) 

A  ( 3 , 3 ) =A3 
A(3,7)  =  -K  1/2.0 
A ( 3 , 11) =-K2/2 . 0 
A ( 4 , 1 )=(P2*S4)/(2.0*M) 

A(4,2)=-P1/2.0 
A ( 4 , 4 ) =A4 
A(4,8)=-K  1/2.0 
A ( 4 , 12 ) =-K2/2 . 0 
A(5, 1 )=-K1/2.0 
A  ( 5 , 5 )  =  N 1 
A( 6 , 2 )  =  -K  1  /2 . 0 
A(  6 , 6 ) =N2 
A! 7 , 3 ) =-K 1 /2 . 0 
A  (  7 , 7  )  =  N  3 
A( 8 , 4 )  =  -K 1 /2 . 0 
A ( 8 , 8 ) =N4 
A ( 9 , 1 )  =  - K2/2 . 0 
A ( 9 , 9 ) =Q1 
At  1 0 , 2 )  =  -K2/2 . 0 
A( 10, 10 ) =02 
A(1 1 , 3 )  =  -K2/2 . 0 
A ( 1 1 , 1 1 ) =Q3 
A ( 12,4) =-K2/2 . 0 
A ( 12, 12) =Q4 
A(  13, 1 3) =B1 
A  (  14 , 14) =B2 
A ( 15, 1 5 ) =C 1 
A ( 1 6 , 16 ) =C2 
A ( 17, 17) =C3 

CALL  EIGRF(A,N, IA, IJOB,W,Z, IZ.WK, IER) 

WRI TE ( 6 , 20 ) 

FORMAT!' 0'  , 12X , ' REAL'  , 1 0X , ' IMAG'  /) 

WRITE ( 6 , 30 )  W 
FORMAT  ( 2X , 2 F 16.4) 

DO  40  1=1,17 
WK ( I ) =REAL ( W ( I  )  ) 

EIGEN  =  AMIN1 ( WK ( 1 ) , WK ( 2 ) , WK ( 3 ) , WK ( 4 ) , WK ( 5 ) , WK ( 6 ) , WK ( 7 ) , WK ( 8 ) , WK ( 9  )  , 
*  WK ( 10) ,WK( 1 1 ) , WK ( 12) , WK ( 13) , WK ( 14) , WK ( 15) , WK ( 16) , WK ( 17) ) 

RETURN 

END 

PRINT  THE  OUTPUT. 

SUBROUTINE  WR 1  TER ( K , A , B , JOPT , TBL ) 

DIMENSION  A ( 4 ) , B ( 4 ) , TBL ( 4 1 , 25 ) 

REAL  K.JOPT 

WRITE ( 6 , 10)  K 

FORMAT!/////  '  K=  ' ,F4.2  /) 


< 
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WR I T  E ( 6 , 15) 

FORMAT!  '  A 1 

WRITE (6, 16)  (A! I ), 1=1 ,4) 

FORMAT! 1H  , 4F 1 2 . 8  /) 


A2 


A3 


A4  '  ) 


WRITE ( 6 , 17) 

FORMAT!  '  B 1 

WRITE ( 6 , 18)  (B(I),I=1,4) 

FORMAT ( 1H  ,  4F 12 . 8  ) 


B2 


B3 


B4  '  ) 


WRITE ( 6 , 19) 
FORMAT!//' 


TIME ( SEC )  Y 1 
Y6 


Y2 


Y3 


Y  4 


Y5 


Y7 


PI 


P2  '  / ) 


T  =  0 . 0 

DO  21  1=1,11 

WRITE ( 6, 20 )  ( T , ( TBL ( I , d ) , J= 1 , 9 ) ) 

FORMAT ( 1 H  , F6 . 2 , 6X , 9F 12 . 8 ) 

T=T+0 . 05 

WRI TE ( 6 , 23 ) 

FORMAT!  /'  TIME! SEC )  P3  P4  P5  P6 

*  P7  P8  P9  P10  P 1 1 ' / ) 

T  =  0 . 0 

DO  25  1=1,11 

WRITE (6, 24)  (T, (TBL! I ,0) , J= 10 , 18)  ) 

FORMAT ( 1 H  ,F6.2,6X,9F 12.8) 

T=T+0 . 05 

WRITE ( 6 , 37  ) 

FORMAT!  /'  TIME! SEC )  P12  P13  P14  P15 

*  VI  V2  EIGENV  '  /) 

T  =  0 . 0 

DO  28  1=1,11 

WRITE (6, 35)  (T,  (TBL! 1 ,0) ,0=19,24)  ) 

FORMAT! 1H  , F6 . 2 , 6X , 6F 1 2 . 8 , F 12 . 4 ) 

T=T+0 . 05 

WR ITE ( 6 , 30 )  OOPT 

FORMAT!  //  '  0 ( CST  )  ='  , F 1 2 . 8  //  ) 

IF!  I ,GT. 1 )  GO  TO  38 
WRITE ( 6 , 19) 

T  =  0 . 5 

DO  32  1=11,41 

WRITE ( 6 , 20 )  (T, (TBL (1,0), 0=1, 9)) 

T=T+0 . 05 

WRITE (6,23) 

T  =  0 . 5 

DO  34  1=11,41 

WRITE ( 6 , 24 )  (T,  (TBL(I,0) ,0=10, 18) ) 

T=T+0 . 05 


WRITE ( 6 , 37 ) 

T  =  0 . 5 

DO  36  1=11 ,41 


' 
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WRITE ( 6 , 35 )  (T,(TBL(I,3),3=19,24)) 

5  FORMAT ( 1 H  , F6 . 2 , 6X , 6F 1 2 . 8 ) 

6  T=T+0 . 05 

7  FORMAT (  /'  TIME ( SEC )  P12  P13  P14 

*  VI  V2  ' / ) 


8  RETURN 
END 

RUNGE-KUTTA  METHOD. 

SUBROUTINE  RUNTA ( N , K , I , X , DX , T , H ) 
DIMENSION  Y( 100) ,Z( 100) , X ( 88 ) , DX ( 88 ) 


K=K+  1 

GO  TO  (1,2,3,4,5),K 

DO  10  3=1 ,N 
Z  (  3 ) =  DX ( J ) 

Y ( 3 ) =X ( 3 ) 

0  X(3)=Y(3)+0.5*H*DX(3) 

5  T=T+0 . 5*H 

1=1 

RETURN 

DO  15  3=1  ,N 
Z(3)=Z(3)+2.0*DX(3) 

5  X(3)=Y(3)+0.5*H*DX(3) 

1  =  1 

RETURN 

DO  20  3=1 ,N 
Z(3)=Z(3)+2.0*DX(3) 

0  X(3)=Y(3)+0.5*H*DX(3) 

GO  TO  25 

DO  30  3=  1  ,N 

0  X(3)=Y(3)  +  (Z(3)+DX(3)  )  *  H  /  6 . 0 


1=2 
K  =  0 

RETURN 

END 

THE  MILNE'S  FORMULA  TO  PREDICT  X(K+1). 

SUBROUTINE  PMI LNE ( M , N , X , DX , T , H ) 

DIMENSION  X ( N ) , DX ( N ) 

DO  1  I  =  1  ,  M 
XT  =  X ( 3*M+I  ) 

X ( 3*M+I ) =X ( 2*M+I ) 

X ( 2*M+ 1 ) =X ( M+I  ) 

X ( M+I ) =X ( I  ) 

X(I)=XT+4.0*H*(2.0*DX(I )-DX(M+I ) +2 . 0*DX ( 2*M+I )  )/3.0 
DX ( 2*M+I ) =DX ( M+I  ) 

DX ( M+I ) =DX ( I ) 


P 1 5 


T  =  T+H 
RETURN 


* 


. 


END 
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THE  MILNE'S  FORMULA  TO  CORRECT  X(K+1) 

SUBROUTINE  CMI LNE ( M , N , X , DX , T , H ) 

DIMENSION  X ( N ) , DX ( N ) 

DO  2  1=1 ,M 

X(I)=X(2*M+I)+H*(DX(I)+4.0*DX(M+I )+DX(2*M+I ) 1/3.0 

RETURN 

END 

THE  EVALUATION  OF  THE  COST  FUNCTIONAL  EQUATION. 

FUNCTION  COST ( Y , A , B , V ) 

DIMENSION  Y(88),A(4),B(4),V(2) 

REAL  N 1 ,N2,N3,N4 

DATA  A1/2.5/.A2/1 .0/.A3/0. 1 / , A4/2 . 5/ , B 1 / 1 .0/.B2/1 .0/.C1/0.01/, 
*  C2/0.01/ ,03/0.01/ 

DATA  N1 ,N2,N3,N4,Q1 , Q2 , Q3 , Q4/8*0 . 00 1 / 


COST  =  A 1 *Y ( 1 )**2+A4*Y(4)**2+A2*Y(2)**2+A3*Y(3)**2+C1*Y( 51**2+ 

*  C2*Y(6)**2+C3*Y(7I**2+B1*V( 1 )**2+B2*V(2)**2+N1*A( 1 1**2  + 

*  N2*A(2)**2+N3*A(3)**2+N4*A(4)**2+Q1*B( 1 ) **2+Q2*B ( 2 1 **2+ 

*  Q3*B(3)**2+Q4*B(4)**2 
RETURN 

END 

INTEGRATION  BY  SIMPSON'S  RULE. 

SUBROUTINE  SIMP ( F , AREA  1 
DIMENSION  F ( 3 1 ) 

REAL  H , XMAX , XMI N 

DATA  H/0 . 05/ , XMAX/0 ,5/,XMIN/0.0/,N/11/ 

H= ( XMAX-XMIN ) / ( N- 1 ) 

X=XMI N+H 
SUM=0 . 0 

DO  4  I = 1 , N 
IF (MODI  1,2)1  2,2,3 
SUM  =  SUM+4 . 0*F (  I  1 
GO  TO  4 

SUM=SUM+2 . 0*F ( 1 1 
X  =  X+H 

AREA= ( H/3 . 0 1  * ( F ( 1 )+SUM+F(N) 1 

RETURN 

END 

TO  SOLVE  VI , V2 , K 1  AND  K2 . 

SUBROUTINE  VVKK ( Y , A , B , V , K 1 , K2 1 
DIMENSION  Y(88) , A ( 4 ) , B ( 4 1 , V ( 2  1 
REAL  K 1 , K2 

DATA  G2/1 . 33/ ,G3/1 .42/.CR/1 . 0/ , TG/0 . 2/ , TE/0 . 2/ , B 1 / 1 .0/ ,B2/1 .0/ 

V( 1 )=A( 1 )*Y( 1  )  +  A(2)*Y(2)+A(3)*Y(3)+A(4)*Y(4) 

V ( 2 ) =B ( 1 )*Y( 1 )+B(2)*Y(2)+B(3)*Y(3)+B(4)*Y(4) 


■ 

9 


K1  =  (Y(22)*G2*G3)/TG-2.0*B1*V(  1  ) 
K2=( Y (21 )*CR)/TE-2.0*B2*V(2) 
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RETURN 

END 


GRAPH  OF  DELTA  VS  TIME  AT  K=1.0 
PROGRAM  PLOT  1 

DIMENSION  X ( 1 1 ) , Y ( 1 1 ) , ALPH ( 20 ) , DATA  1 ( 11,2) 

REAL  ALPH/'  FIG.'  '  ,'  ROTO'  ,'  R-AN'  ,'  GLE  '  , '  VS  T'  ,'  IME  '  ,'  (  AT'  , 

*/  K=  '  '  1  0  '  '  )  '  '  ' 

*'  TIME'  !'  ( SE'  !'  COND'  j'  )  '  /  ROTO'  ,'  R-AN'  , '  GLE  ( '  ,'  DEG)'  / 

DATA  DATA  1/37. 1,38.4,42.1,48.0,55.2,62.7,70.0,75.9,79.7,81.5, 

*81 . 1 ,37. 1 ,38.4,42. 1 ,47.9,54.8,61 .9,68.0,72.0,73.5,72.4,68.6/ 

ND=  1 1 
NF=  1 
KC  =  3 
HA  =  0 . 0 
HB  =  0 .  1 
HC=5 . 0 
VA=0 . 0 
VB=  1  0 
VC=8 . 5 


X  (  1 )  =0 . 0 

DO  10  1=2,11 

X  ( I )  =X ( I -  1 ) +0 . 05 

DO  30  1=1,2 
DO  20  J=  1  ,  1  1 
Y  (  J ) =DATA 1 ( J , I ) 

CALL  CGPL2(X, Y ,ND,NF , KC , HA , HB , HC , VA, VB, VC, ALPH) 

NF  =  -2 

CONTINUE 

STOP 

END 


GRAPH  OF  DELTA  VS  TIME  AT  K=0.4 
PROGRAM  PLOT2 

DIMENSION  X (  1 1  )  ,  Y ( 1 1 ) , ALPH ( 20 ) , DATA2 ( 11,2) 

REAL  ALPH/'  FIG.'  ,'  '  ,'  ROTO'  ,'  R-AN'  ,'  GLE  '  , '  VS  T'  ,'  IME  '  ,'  (  AT'  , 

*'  K=  '  '  0 . 4  '  '  )  '  '  ' 

*'  TIME'  /  ( SE'  !'  COND'  \  ‘  ROTO'  ,'  R-AN'  , '  GLE  ( '  ,'  DEG)'  / 

DATA  DATA2/37. 1,37.6,39.1,41.5,44.3,47.5,50.5,53.1,55.0,56.0, 

*  56.1,37.1,37.6,39.1,41.5,44.2,47.1,49.5,51.3,52.1,51.7,50.3/ 


ND=  1 1 
NF=  1 
KC  =  3 
HA  =  0 . 0 
HB  =  0 .  1 
HC=5 . 0 
VA=0 . 0 
VB=  1  0 
VC=8 . 0 


m  in  *  r;  :: 
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X(1)=0.0 
DO  10  1=2,11 
0  X( I )=X(  1-1 1+0.05 

DO  30  1=1,2 
DO  20  0=1,11 
0  Y ( J ) =DATA2 (0,1) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC, ALPH) 
NF  =  -2 

D  CONTINUE 
STOP 
END 


GRAPH  OF  DELTA  VS  TIME  AT  UNSTABLE  POINT 
PROGRAM  PLOT 3 

DIMENSION  X ( 1 1 ) , Y ( 1 1 ) , ALPH (20) , DATA3 ( 11,2) 

REAL  ALPH/'  FIG.'  ,'  '  ROTO'  , '  R-AN'  ,'  GLE  '  VS  T'  ,'  IME  '  ,'  (  AT'  , 

*'  K=  '  '  2  1  '  1  )  11  1 

♦  'TIME'  ! '  isE'  !' COND'  !'  )  '/  ROTO'  ,'  R-AN'  ,'  GLE  (','  DEG )'  / 

DATA  DATA3/37 . 1,39.9,47.8,60.0,75.0,90.0,104.3,117.1,128.0, 
*137.5,146.9,37.1,39.9,47.8,59.8,73.9,90.0,98.9,107.2,110.8, 

*109.8, 104.3/ 


ND=  1 1 
NF=  1 
KC=3 
HA  =  0 . 0 
HB  =  0 . 1 
HC  =  5 . 0 
VA  =  0 . 0 
VB  =  20 
VC=8 . 0 


X(  1 ) =0 . 0 
DO  10  1=2,11 
3  X( I ) =X ( I  -  1 )  +  0 . 05 

DO  30  1=1,2 
DO  20  J=  1 ,  1 1 
3  Y(J)=DATA3(d,  I) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA, VB,VC, ALPH) 
NF  =  -2 

3  CONTINUE 
STOP 
END 


GRAPH  OF  DELTA  VS  TIME  FOR  K=0.5  BUT  USE  OPTIMAL  AI.BI  FOUND 
FOR  K=0 . 2 


PROGRAM  PLOT 4 

DIMENSION  X ( 1 1 ) , Y ( 1 1 ) , ALPH ( 20 ) , DATA4 ( 1 1,2) 

REAL  ALPH/'  FIG.'  ,'  '  ,'  ROTO'  ,'  R-AN'  ,' GLE  '  ,'  VS  T'  ,'  IME  '  ,'  ( 

*'  K=  '  , '  0 . 5  '  , '  )  '  '  ' 

♦'TIME'  ] '  ( SE'  !'  COND'  /  )  '/ ROTO'  ,' R-AN'  ,' GLE  (',' DEG)' / 

DATA  DATA4/37. 1 ,37.7,39.6,42.5,46.0,49.5,52.6,54.9,55.8,55.3, 
*  53.4,37.1,37.7,39.6,42.5,46.0,49.5,52.6,54.7,55.5,54.5, 


AT' 


8H  AH  n .  .  OH , 


52.7/ 
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ND=  1 1 
NF=  1 
KC=3 
HA=0 . 0 
HB=0 . 1 
HC=5 . 0 
VA=0 . 0 
VB=  1 0 
VC-8 . 0 


X  ( 1 )  =0 . 0 
DO  10  1=2,11 
10  X  ( I )  =X  ( I  - 1  )  +0 . 05 

DO  30  1=1,2 
DO  20  0=1,11 
>0  Y(  d )  =DATA4  (  d  ,  I  ) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC,ALPH) 
NF  =  -2 

10  CONTINUE 
STOP 
END 


GRAPH  OF  DELTA  VS  TIME  FOR  K=0.1  T=0.0  TO  T=2.0 
PROGRAM  PLOT5 

DIMENSION  X ( 2 1 ) , Y ( 2 1 ) , ALPH ( 20 ) , DATA5 ( 21 ,2) 

REAL  ALPH/'  FIG. '  ROTO'  ,'  R-AN'  ,'  GLE  '  VS  T'  ,'  I  ME  '  , ' 

*'  K=  ','0.1  '  , '  )  '  , '  '  , 

*'  TIME'  , '  { SE'  ,'  COND'  '/  )  ',' ROTO'  ,' R-AN'  ,' GLE  (',' DEG)' / 

DATA  DAT A5/37 . 1,37.6,38.9,40.3,41.2,41.1,39.9,37.7,35.6,34.2 
*34.0,35.0,36.8,38.5,39.6,39.7,38.8,37.4,36.0,35.1,35.2, 
*37.1,37.6,38.9,40.5,41.7,42.1,41.2,39.0,36.3,34.2,33.5, 
*34.3,36.0,38.2,39.8,40.4,39.7,38.2,36.4,35.2,34.7/ 


ND=2 1 
NF=  1 
KC=3 
HA  =  0 . 0 
HB  =  0 . 4 
HC=5 . 0 
VA=0 . 0 
VB=  1 0 
VC=7 . 0 


X( 1 ) =0 . 0 
DO  10  1=2,21 
1°  X(  I )  =X  ( I  -  1  )  +0 .  1 

J 

DO  30  d=1 ,2 
DO  20  1=1,21 
Y  ( I )  =DAT  A5 1 1 ,  d ) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,V6,VC,ALPH) 

NF  =  -  2 

30  CONTINUE 
STOP 
END 


(  AT'  , 
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GRAPH  OF  PARAMETERS! AI )  VS  DISTURBANCE 
PROGRAM  PL0T6 

DIMENSION  X (  1 1  )  ,  Y ( 1 1 ) , ALPH ( 20  )  , DATA6 ( 1 1,4) 

REAL  ALPH/'  FIG.'  '  , '  OPT  I '  , '  MAL  '  PARA'  , '  METE'  ,'  RS  V'  , 

*'  S  DI '  ,'  STUR'  ,'  BANC'  ,'  E  '  DIST'  , '  URBA'  , '  NCE  '  , 

*'  (K)  '  ,'  PARA'  ,'  METE'  RS  '  ,'  (Alt'  / 

DATA  DATA6/0. 0,-0. 01 2, 0.0 16 ,0.094,0. 172,0.196,0.224,0.250,0.273, 
*0.290,0.2999,0.0,-0.242,-0.419,-0.420,-0.416,-0.407,-0.399,-0.392, 
*-0.387,-0.382,-0.378,0.0,0.00,-0.001,-0.002,-0.003,-0.002,-0.001, 
*0.001,0.003,0.007,0.01,0.0,0.009,-0.016,-0.086,-0.16,-0.184, 
*-0.213,-0.242,-0.268,-0.288,-0.303/ 


ND=  1 1 
NF  =  1 
KC=3 
HA=0 . 0 
HB  =  0 . 2 
HC  =  5 . 0 
VA=-0 . 5 
VB  =  0 .  1 
VC=8 . 0 


X  ( 1 )  =0 . 0 
DO  10  1=2,11 
X ( I )=X( 1-1 )+0. 1 

DO  30  1=1,4 
DO  20  J=1  ,  1 1 
Y I J ) =DATA6 ( J  ,  I  ) 

NF  =  I 

I F ( NF . EQ . 4 )  NF  =  -4 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC,ALPH) 

CONTINUE 

STOP 

END 


GRAPH  OF  PARAMETER ( B I )  VS  DISTURBANCE 
PROGRAM  PL0T7 

DIMENSION  X< 1  1  )  ,  Y( 1 1 ) , ALPH (20) ,DATA7 ( 11,4) 

REAL  ALPH/'  FIG.'  ,'  '  , '  OPT  I'  , '  MAL  '  , '  PARA'  ,'  METE'  ,'  RS  V'  , 

*'  S  DI'  ,'  STUR'  ,'  BANC'  , '  E  '  ,'  '  ,'  DIST'  ,'  URBA'  , '  NCE  '  ,'  (K)  ' 

*'  PARA'  ,'  METE'  ,'  RS  '  ,'  (  B I  )'  / 

DATA  DATA7/0 . 0  ,  0.0032,-0.0045,-0.0258,-0.048,-0.055,-0.0637, 
*-0.0724, 

*-0.0804,-0.0867,-0.0917,0.0,0.065,0. 1138,0.1163,0.1181,0.1174, 
*0.1176,0.1180,0.1185,0.119,0.1198,0.0,0.00,0.0,0.0,0.0,0.0, 
*-0.0003,-0.0007,-0.0015,-0.0023,-0.0032,0.0,-0.0024,0.0043, 
*0.0235,0.0446,0.0516,0.0606,0.0702,0.0796,0.0872,0.0936/ 


ND=  1 1 
N  F  =  1 
KC=3 
HA=0 . 0 
HB  =  0 . 2 
HC  =  5 . 0 


VA=-0 .  12 
VB=0 . 03 
VC=8 . 0 
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X( 1  )  =  0. 0 
DO  10  1=2,11 
10  X(  I )  =X  ( I  -  1  )  +0 .  1 

DO  30  1=1,4 
DO  20  0=1,11 
>0  Y(J)=DATA7(0,I  ) 

NF  =  I 

IF ( NF . EQ . 4 )  NF  =  -4 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC,ALPH) 
!0  CONTINUE 
STOP 
END 


GRAPH  OF  OPTIMUM  CONTROLS  VS  TIME  AT  K=0.4 
PROGRAM  PL0T8 

DIMENSION  X(11),Y(11), ALPHA (20) ,ALPHB(20) , DATA8 ( 11,2) 

REAL  ALPHA/'  FIG.'  ,'  '  , '  OPT  I '  , '  MUM  '  , '  CONT'  ,'  ROL  '  ,'  TO  T'  , 

*'URBO'  ,'  -GEN'  ,'  ERAT'  , '  OR  M'  ,'ODEL'  ,'  TIME'  ,'  (  SE'  , '  COND'  , '  ) 

*'  SPEE'  ,'  DER  '  ,'  GEAR'  , '  '  / 

REAL  ALPHB/'  FIG.'  ,'  '  , '  OPT  I '  , '  MUM  '  , '  CONT'  ,'  ROL  '  ,'  TO  T'  ,'  URBO' 

*,' -GEN'  ,' ERAT'  ,' OR  IT  ,'ODEL'  ,' TIME'  ,'  (  SE'  , '  COND'  , '  ) 

* /  E  XC 1 7  '  TER  '  '  F  I  EL'  'D  '  / 

DATA  DATA8/0. 424, 0.273,0. 148,0.065,0.035,0.066,0.156,0.293,0.456, 
*0.621,0.759,0.738,0.780,0.816,0.840,0.848,0.840,0.814,0.775,0.730, 
*0.683,0.644/ 


ND=  1 1 
NF=  1 
KC  =  3 
HA  =  0 . 0 
HB=0 . 1 
HC=5 . 0 
VA=0 . 0 
VB  =  0 . 2 
VC=4 . 0 


X  ( 1 )  =0 . 0 
DO  10  1=2,11 
0  X ( I ) =X ( I  -  1 ) +0 . 05 

DO  20  1=1,11 
0  Y ( I ) =DATA8 (1,1) 

CALL  CGPL2(X, Y,ND,NF , KC , HA , HB , HC , VA , VB , VC, ALPHA) 

NF  =  - 1 
VA=0 . 6 
VB=0 . 1 
VC=3 . 0 

A  DO  30  1=1,11 
0  Y ( I ) =DATA8 (1,2) 


0,£V  I.O.t  Jft.C  t'^l 


CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB, VC, ALPHB) 

STOP 

END 
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•/ 


C 


10 


20 


r 


GRAPH  OF  MECHANICAL  TORQUE  AND  FIELD  VOLTAGE  DUE  TO  OPTIMUM 
CONTROLS  AT  K  =  0.4 
PROGRAM  PLOT9 

DIMENSION  X( 1 1  )  , Y ( 1 1 ) , ALPHA <  20 ) , ALPHB ( 20 ) , DATA9 ( 1 1,2) 

REAL  ALPHA/'  FIG.'  '  , '  VARI'  , '  ATIO'  , '  N  OF'  , '  OPT', 

*'  IMUM'  , '  MEC'  ,'HANI'  ,'CAL  '  ,'TORQ'  ,'UE  ','TIME','  (  SE'  , 

*'  COND'  , '  )  '  ,'MECH'  ,' ANIC'  ,' AL  T'.'ORQU'/ 

REAL  ALPHB/'  FIG.'  ,'  ',' VARI'  ,' ATIO'  ,' N  OF'  , '  OPT', 

*'  IMUM'  , '  FIE'  ,' LD  V'  ,'OLTA'  ,' GE  ','  ','TIME','  (SE', 

*'  COND'  , '  )  '  FIEL'  ,'D  VO'  ,' LTAG'  ,'E  '/ 

DATA  DATA9/0. 8, 0.799, 0.792, 0.777, 0.752, 0.720, 0.684, 0.651  , 
*0.626,0.614,0.617,0.738,0.743,0.755,0.771,0.788,0.801, 
*0.807,0.804,0.792,0.773,0.748/ 


ND  =  1 1 
NF=  1 
KC  =  3 
HA=0 . 0 
HB=0.  1 
HC=5 . 0 
VA=0 . 6 
VB=0 . 05 
VC  =  4 . 0 


X ( 1 ) =0 . 0 

DO  10  1=2,11 

X(  I )  =X ( I  - 1 ) +0 . 05 

DO  20  1=1,11 
Y  ( I )  =DATA9 (1,1) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB, VC, ALPHA) 


NF  =  -  1 
VA=0 . 70 
VB=0 . 05 
VC  =  3 . 0 


DO  30  1=1,11 
Y ( I ) =DATA9 (1,2) 


CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB, VC, ALPHB) 

STOP 

END 

r 

is 


GRAPH  OF  OPTIMUM  ANGULAR  VELOCITY  AT  K=0.4 
PROGRAM  PLOT  1 0 

DIMENSION  X( 1 1 ) , Y(  1 1  )  ,  ALPH ( 20 )  ,DATA10(  1  1  ) 

REAL  ALPH/' FIG.'  ,'  ',' VARI'  ,' ATIO'  ,' N  OF'  ,'  OPT', 'IMUM', 

*'  ANG'  ,'ULAR'  ,'  VEL'  ,'OCIT'  ,'Y  ','TIME','  (  SE'  , '  COND'  , 

*'  >  '  ,'ANGU'  ,' LAR  '  ,'  VELO'  , '  CITY'  / 

DATA  DATA  10/0. 0,0. 366, 0.682, 0.907, 1 .006,0.960,0.768,0.456, 


■ 


*0.069,-0.329,-0.675/  143 

C 

ND=  1 1 
NF  =  - 1 
KC = 3 
HA  =  0 . 0 
HB  =  0 .  1 
HC  =  5 . 0 
VA=-0 . 75 
VB=0 . 25 
VC=8 . 0 

C  X( 1 ) =0 . 0 

DO  10  1=2,11 
10  X ( I ) =X ( I  -  1 )+0.05 
C 

DO  20  1=1,11 
20  Y ( I ) =  D A  T  A  10(1) 

C 

CALL  CGPL2 ( X , Y , ND , NF , KC , HA , HB , HC , VA, VB, VC, ALPH) 

STOP 

END 

C 

c 

c 

C  GRAPH  OF  OPTIMUM  EXCITATION  VOLTAGE  AT  K=0.4 
C  PROGRAM  PLOT  11 

DIMENSION  X ( 1 1 ) , Y ( 1 1 ) , ALPH ( 20 ) , DATA  1 1(11) 

REAL  ALPH/'  FIG.7  ,7  7  , ' VARI7  , 7  ATIO7  , 7  N  OF'  , '  OPT7  ,7  IMUM7  , 

*'  EXC7  ,7  I  TAT7  , 7  ION  7  , 7  VOLT7  , 7  AGE  7  ,7  TIME7  , 7  (SE7  , 7  COND7  , 

*7  )  7  ,7  EXCT  ,7  TATI7  ,7  ON  7  , 7  EMF  7  / 

DATA  DATA  1 1/2.889,2.908,2.955,3.018,3.085,3. 135,3.159, 

*3.147,3. 100,3.026,2.928/ 

C 

ND=  1 1 
NF  =  - 1 
KC  =  3 
HA=0 . 0 
HB  =  0 .  1 
HC=5 . 0 
VA=2 . 8 
VB=0 . 05 
VC=8 . 0 
C 

X( 1 )=0.0 
DO  10  1=2,1  1 
10  X ( I ) =X ( I  -  1 ) +0 . 05 
DO  20  1=1,11 
20  Y ( I ) =DAT A 1  1(1) 

C 

CALL  CGPL2(X, Y , ND , NF , KC , HA , HB , HC , VA , VB, VC, ALPH) 

STOP 

END 

C 

C 

c 

S'  GRAPH  OF  DELTA  VS  TIME  FOR  K  =  0 . 4  T  =  0.0  T  =  2.0 
L  PROGRAM  PLOT  1 2 

DIMENSION  X ( 2 1 ) , Y ( 2 1 ) , ALPH (20) , DATA  1 2 ( 2 1 , 2 ) 

REAL  ALPH/7  FIG.7  ,7  7  ,7  ROTO7  ,7  R-AN7  ,7  GLE  7  ,7  VS  T7  ,7  IME  7  , 


' 


C_>  CM  CO  uuuuo 


/ 


/ 


/ 
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*'  (  AT  , '  K=  '  , '  0 . 4  '  , '  ) 

•  'TIME'!'  (SE'  !'  COND'  !'  )  '!'  ROTO'  R-AN'  ,'  GLE  (','  DEG )'  / 

DATA  DAT A 12/37. 1 ,39.1,44.2,49.5,52.1,50.3.44.3,36.1,29.2,25.6, 
•26.4,30.9,37.3,43.2,46.5,45.7,41.8,36.6,32.0,29.8,30.4, 
*37.1,39.1,44.3,50.5,55.0,56.1,51.6,42.3,32.3,25.2,23.0, 
*25.9,32.8,41.0,47.4,49.5,46.8,40.9,34.4,29.7,28.2/ 

C 

ND=2 1 
NF=  1 
KC=3 
HA=0 . 0 
HB=0 . 4 
HC=5 . 0 
VA=0 . 0 
VB  =  1 0 . 0 
VC  =  7 . 0 

c 

X  ( 1 ) =0 . 0 
DO  10  1=2,21 
10  X ( I ) =X ( I  -  1 ) +0 .  1 

DO  30  1=1,2 
DO  20  J= 1 ,21 
Y(  d ) =DAT A 1 2 ( J , I ) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC,ALPH) 

NF  =  -2 
CONTINUE 
STOP 
END 


C 


c 


PHASE  PLANE  PLOT  AT  K  =  2. 1 
PROGRAM  PLOT  1 3 

DIMENSION  X  ( 2 1  )  ,  Y  (  2 1 ) , AL  PH ( 20 ) , DATA  13(21 ,4) 

REAL  ALPH/' FIG.'  ,'  '.'PHAS'.'E  PL'.'ANE  ','PLOT', 

*'  AT  '  , '  K=  ','0.4  '  , '  '  , '  '  , '  '  , 

*'  ROTO'  !'  R-AN'  ! '  GLE  ( '  !'  DEG)'  !'  ANGU'  !'  LAR-'  /  VELO'  , '  CITY'  / 

DATA  DATA  13/37. 1 ,39. 1,44.2,49.5,52.1,50.3,44.3,36.1,29.2, 
*25.6,26.4,30.9,37.3,43.2,46.5,45.7,41.8,36.6,32.0,29.8,30.4, 
*37.1,39.1,44.3,50.5,55.0,56.1,51.6,42.3,32.3,25.2,23.0, 
*25.9,32.8,41.0,47.4,49.5,46.8,40.9,34.4,29.7,28.2, 
*0.0,0.68,1.01,0.77,0.07,-0.67,-1.34,-1.40,-0.96,-0.24,0.49, 
*1.02,1.15,0.85,0.22,-0.45,-0.87,-0.91,-0.62,-0.14,0.34, 
*0.0,0.69,1.07,0.99,0.51,-0.14,-1.32,-1.82,-1.57,-0.83,0.09, 
*0.91,1.39,1.36,0.79,-0.07,-0.82,-1.17,-1.03,-0.56,0.06/ 


ND  =  2 1 
NF  =  1 
KC=3 
HA= 10.0 
HB= 10.0 
HC  =  5 . 0 
VA  =  -2.0 
VB  =  0 . 5 
VC=8 . 0 


10 


DO  30  d=  1 , 2 
DO  10  1=1,21 
MI)  =DATA  1 3  ( I ,  J  ) 


, 
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DO  20  1=1,21 
0  Y ( 1 1 =DATA 1 3(1, d  +  2 ) 

CALL  CGPL2(X,Y,ND,NF,KC,HA,HB,HC,VA,VB,VC,ALPH) 
NF  =  -2 

0  CONTINUE 
STOP 
END 


PHASE  PLANE  PLOT  AT  K  =  2. 1 
PROGRAM  PLOT  1 4 

DIMENSION  X ( 1 1 ) , Y  < 1 1 ) , ALPH ( 20 ) , DATA  1 4 ( 11,4) 

REAL  ALPH/'  FIG.'  '  ,'  PHAS'  , '  E  PL'  ,'  ANE  '  ,'  PLOT'  ,'  AT  '  , 

K=  '  1  2  1  '  '  11  11  ' 

*'  ROTO'  /  R:AN'  /  GLE('  !'  DEG)'  /  ANGU'  /  LAR-'  , '  VELO'  CITY'  / 

DATA  DATA  14/37. 1,39.9,47.8,59.8,73.9,90.0,98.9, 

*107.2, 110.8,109.8,104.3, 

*37.1,39.9,47.8,60.0,75.0,90.0,104.3,117.1,128.0,137.5,146.9, 

*0.0,1.92,3.57,4.68,5.04,4.59,3.49,2.03,0.43,-1.14,-2.55, 

*0.0,1.93,3.60,4.81,5.40,5.32,4.77,4.07,3.49,3.22,3.39/ 


ND=  1 1 
NF=  1 
KC  =  3 
HA  =  0 . 0 
HB  =  30 . 0 
HC=5 . 0 
V A=-3.0 
V  B  =  1  .0 
VC=8 . 5 

DO  30  J=1  ,2 
DO  10  1=1,11 
0  X ( I ) =DATA 14(1, J) 

DO  20  1=1,11 
0  Y ( I ) =DATA 14(1, J  +  2  ) 

CALL  CGPL2(X, Y,ND,NF , KC , HA , HB , HC , VA, VB, VC, ALPH) 
NF  =  -2 

0  CONTINUE 
STOP 
END 


* 


■ 

•- 


